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SUMMARY 


This  report  summarizes  the  results  of  the  research  in  two  parts.  Part  I  describes  the 
development  and  application  of  a  free-wake  Euler  and  Navier-Stokes  computational  fluid 
dynamics  method,  called  ’’TURNS”,  for  helicopter  applications.  This  finite-difference, 
implicit,  upwind  numerical  method  uses  structured  grids,  and  has  been  used  for  calculating 
the  viscous,  three-dimensional  flowfields  of  rotors  in  hover,  forward  flight,  blade-vortex 
interactions,  and  high  speed  impulsive  noise.  In  this  free-wake  method,  the  induced  effects 
of  wake,  including  the  interaction  of  tip  vortices  with  successive  blades,  are  captured  as  a 
part  of  the  overall  flowfield  solution  without  specifying  any  wake  models. 

The  hovering  flowfield  calculations  are  done  for  a  two-bladed  model  rotor  as  well  as 
for  four-bladed  modern  rotors  of  UH-60A  Blackhawk  and  BERP  helicopters.  Different 
modifications  of  the  UH-60A  blade  have  also  been  considered.  Using  the  computational 
results,  an  attempt  is  made  to  understand  the  importance  of  planform  effects  of  the  var¬ 
ious  four-bladed  rotor  configurations  considered.  Calculated  results  are  presented  in  the 
form  of  surface  pressures,  hover  performance  parameters,  surface  partice  flow,  tip  vortex 
formation,  and  vortex  wake  trajectory.  Comparison  of  calculated  surface  pressures  and  the 
nearfield  wake  trajectory  for  the  model  two-bladed  rotor  and  the  UH-60A  Blackhawk  ro¬ 
tor  show  good  agreement  with  measured  data.  The  captured  vortex  structure  is,  however, 
diffused  due  to  coarse  grids,  but  this  appears  to  have  minimal  influence  on  the  prediction 
of  surface  pressures  and  airloads.  Comparison  of  UH-60A  results  with  an  equivalent  rect¬ 
angular  UH-60A  blade  and  a  high-twist  BERP  blade  indicates  that  the  BERP  blade,  with 
an  unconventional  planform,  produces  more  thrust  at  a  given  collective  pitch,  and  approx¬ 
imately  the  same  Figure  of  Merit.  The  high  thrust  conditions  considered  produce  severe 
shock-induced  flow  separation  for  the  UH-60A  blade,  while  the  BERP  blade  develops  more 
thrust  and  minimal  separation.  The  BERP  blade  produces  a  tighter  tip  vortex  structure 
compared  with  the  UH-60A  blade.  These  results  and  the  discussion  presented  bring  out 
the  similarities  and  differences  between  the  various  four-bladed,  high  performance  rotors 
considered. 
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The  capabilities  of  the  TURNS  code  are  further  demonstrated  by  calculating  the  blade- 
vortex  interactional  flowfield  of  a  rotor  in  forward  flight  and  encountering  an  upstream 
generated  concentrated  line  vortex.  The  code  has  also  been  used  to  calculate  the  acoustics 
of  high  speed  impulsive  noise  in  both  hover  and  forward  flight.  Sample  results  are  presented 
for  these  cases  and  compared  with  available  experimental  data.  Limited  comparisons  of 
the  Navier- Stokes  results  with  the  Euler  results  are  also  presented  for  both  aerodynamics 
and  acoustics. 

In  Part  II,  the  unsteady  flowfield  results  of  an  oscillating  wing  are  presented.  The 
unsteady  flowfields  of  a  two-dimensional  oscillating  wing  are  calculated  using  an  implicit, 
finite-difference,  Navier-Stokes  numerical  scheme  using  five  widely  used  turbulence  models. 
The  objective  of  this  study  is  to  identify  an  appropriate  turbulence  model  for  accurate  sim¬ 
ulation  of  three-dimensional  dynamic  stall.  Three  unsteady  flow  conditions  corresponding 
to  attached  flow,  light-stall,  and  deep-stall  of  an  oscillating  wing  experiment  were  chosen 
as  test  cases  for  computations.  Results  of  unsteady  airload  hysteresis  curves,  harmonics  of 
unsteady  pressures,  and  instantaneous  flowfield  pictures  are  presented. 

Unsteady  airloads  calculated  using  the  turbulence  models  of  the  Baldwin-Lomax, 
Renormalization  Group  Theory,  the  Johnson-King,  the  Baldwin-Barth,  and  the  Spalart- 
Allmaras  are  compared  with  experimental  data  of  Piziali  of  U.  S.  Army  Laboratory.  The 
comparisons  of  unsteady  airloads  show  that  all  models  are  deficient  in  some  sense  and  not 
a  single  model  predicts  all  airloads  consistently  and  in  agreement  with  experiment  for  all 
three  regimes  of  flow  conditions  considered.  Considering  the  overall  performance  of  these 
models,  the  Baldin-Barth  one-equation  model  appears  to  have  a  better  performance  when 
the  flow  is  separated.  It  is  also  found  to  be  the  most  expensive  model  in  terms  of  compu¬ 
tational  cost.  Finally,  the  three-dimensional  oscillating  wing  calculations  done  using  the 
Baldwin-Lomax  turbulence  model  and  a  modified  TURNS  algorithm  for  the  deep-stall  case 
demonstrate  that  the  simple  algebraic  model  consistently  underpredicts  separation  history 
just  like  in  the  two-dimensional  case  and  the  predicted  airloads  have  poor  agreement  with 
experiments. 
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PART  I 

TURNS  CODE  AND  ITS  APPLICATION  TO  HELICOPTER  ROTORS 
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1.  INTRODUCTION 


The  accurate  numerical  simulation  of  the  helicopter  rotor  flowfield  in  hover  continues  to 
be  one  of  the  most  complex  and  challenging  problems  of  applied  aerodynamics.  This  is 
true  despite  the  availability  of  the  present  day  supercomputers  and  improved  numerical 
algorithms.  The  complexity  of  the  flow  stems  from  several  peculiar  features  that  are  unique 
to  the  helicopter  rotor;  the  dominant  of  these  is  the  vortical  wake,  since  is  is  responsible 
for  producing  unsteady  loads  fluctuation,  noise,  and  vibration.  In  other  words,  the  vortical 
wake  strongly  influences  the  operating  characteristics  of  helicopter.  Accurate  prediction 
of  this  vortical  wake  is  probably  the  most  important,  most  studied  and  the  most  difficult 
aspect  of  helicopter  flowfield. 

Current  methods  of  analysis  of  wakes  range  in  complexity  from  relatively  simple 
momentum-theory  applications  to  sophisticated  free-wake  methods.  In  between  these  ex¬ 
tremes,  there  exists  a  variety  of  so-called  prescribed- wake  models1 .  Although  such  methods 
are  often  used  in  industry  for  predicting  airloads,  they  do  suffer  from  the  limitation  that 
the  empirical  determination  of  wake  shape  ignores  some  of  the  important  flowfield  details 
such  as  the  mutual  interaction  of  individual  vortex  elements.  Since  most  of  these  flowfield 
prediction  methods  have  to  be  coupled  with  some  wake  model  for  realistic  estimates  of  air¬ 
loads,  they  are,  therefore,  restricted  in  their  application  to  blade  shapes  and  planforms  that 
are  geometrically  simple.  Besides,  further  modeling  is  required  to  simulate  the  formation 
of  concentrated  tip  vortex  in  the  tip  region. 

Although,  a  complete  flowfield  simulation  of  helicopter  is  still  beyond  the  current  ca¬ 
pability,  a  combination  of  powerful  supercomputers  in  conjuction  with  improved  numerical 
algorithms  have  enabled  advances  to  be  made,  recently,  using  Computational  Fluid  Dynam¬ 
ics  (CFD)  to  solve  the  equations  of  fluid  motion  for  individual  elements  of  such  a  complex 
flowfield.  The  solution  schemes  for  these  equations  are  usually  coupled  with  integral  wake 
models,  e.  g.  vortex  line  elements,  vortex  lattices,  or  panels,  to  bring  in  the  influence  of 
the  vortex  wake.  Techniques  that  implement  this  idea  using  a  “prescribed”  wake  geometry 
encompass  potential  flow1-5,  Euler6-8  and  Navier-Stokes  methods9-12.  Also,  a  potential 
flow  method13  has  been  coupled  with  a  “free”  wake  approach,  in  which  the  wake  vorticity  is 
specified,  but  is  allowed  to  convect  freely  with  the  flow  without  constraining  its  trajectory. 
Other  recent  methods,  including  the  present  one,  use  direct  techniques  for  the  entire  solu¬ 
tion  process,  using  Euler14-17  and  Navier-Stokes18  flow  solvers  by  capturing  the  influence 
of  wake  as  a  part  of  the  overall  flowfield  solution.  In  this  context  one  can  refer  to  these 
“free-wake”  methods  as  wake-capturing  schemes,  in  analogy  to  shock-capturing,  whereas 
the  “coupled-”  or  “prescribed-wake”  methods  are  somewhat  analogous  to  shock-fitting.  In 
another  study,  a  totally  different  approach  has  been  followed  that  is  different  from  the 
above  schemes,  viz.,  the  Navier-Stokes  calculations  have  been  performed19  for  the  global 
time-averaged  wake  alone,  with  the  details  of  the  flow  on  the  rotor  itself  prescribed. 

The  weakest  link  in  the  wake-coupled  methodologies1-12  is  the  wake  modeling.  That 
is,  the  technique  of  specifying  a  prescribed  wake  has  to  be  specialized  for  each  blade  shape, 
making  it  difficult  to  treat  blade  shapes  with  arbitrary  twist,  taper,  and  planforms.  There- 


2 


fore  ,  the  purpose  of  this  study  is  to  develop  an  improved  calculation  method  for  the 
solution  of  Navier-Stokes  equations  for  the  complete  flowfield  of  a  lifting  rotor,  including 
the  wake  and  its  induced  effects.  The  vortex  wake  and  its  effects  are  captured  as  a  part 
of  the  complete  flowfield,  and  thus  no  arbitrary  inputs  or  vortex- modeling  approximations 
are  necessary  to  describe  the  wake.  In  addition  to  the  wake-capturing  capabilities,  the 
Navier-Stokes  approach  is  chosen  for  the  following  reasons:  1)  better  tip-fiow  simulation, 
which  involves  resolving  the  flow  separation  and  the  formation  of  a  concentrated  tip  vortex, 
2)  accurate  simulation  of  strong  viscous-inviscid  interaction  involving  shock  induced  sepa¬ 
ration  at  high  blade  tip  speeds  and  high  collective  pitch  conditions,  and  3)  future  modeling 
of  retreating  blade  and  dynamic  stall  regimes  in  forward  flight. 

A  new  free-wake  Euler  and  Navier-Stokes  CFD  method,  called  TURNS  (transonic  un¬ 
steady  rotor  Navier-Stokes)  Code,  is  developed  for  helicoter  rotors  in  this  study.  This  is  an 
improvement  of  the  version  that  was  developed  previously  by  the  present  author  in  related 
studies  with  wake  modeling9.  One  fundamental  difference  of  the  present  numerical  scheme 
is  the  use  of  Roe’s  upwinding  for  all  three  coordinate  directions20.  This  feature,  coupled 
with  an  implicit  iterative  procedure  has  produced  a  fast,  efficient,  and  accurate  numerical 
scheme.  In  addition,  a  periodic  boundary  condition17  has  been  implemented  in  the  az¬ 
imuthal  direction,  as  described  later  in  the  text.  Also,  a  new  procedure  for  implementing 
farfield  boundary  conditions  for  a  hovering  rotor  has  been  used.  These  improvements  allow 
the  near  wake  to  be  computed  well  enough  to  approximately  simulate  the  correct  inflow 
through  the  rotor,  thus  obviating  the  need  for  the  ad  hoc  wake  modeling  used  previously9. 
These  additional  changes  in  the  Euler/Navier-Stokes  algorithm  are  based  on  some  of  the 
numerical  procedures  described  in  Ref.  21. 

Radiated  noise  can  severely  restrict  rotorcraft  usage  in  both  civilian  and  military  op¬ 
erations.  Impulsive  noise,  the  sum  total  of  the  blade-vortex  interaction  (BVI)  noise  and 
the  high-speed  impulsive  (HSI)  noise,  forms  the  most  important  component  of  the  radiated 
noise.  The  BVI  noise  is  generated  due  to  the  interaction  of  the  vortical  wake  with  rotating 
blades  and  is  generally  more  difficult  to  model  due  to  the  importance  of  unsteady,  three- 
dimensional,  and  vortex  wake  effects.  High-speed  impulsive  noise,  on  the  other  hand,  is 
caused  primarily  due  to  compressibility  effects.  If  the  advancing  blade  tip  Mach  number 
is  highly  supercritical,  the  phenomenon  known  as  delocalization  may  occur,  wherein  the 
supersonic  pocket  on  the  rotor  blade  extends  out  to  the  far-field  beyond  the  blade  tip.  If 
this  occurs,  the  noise  becomes  more  impulsive  and  in  particular  gets  focused  in  its  direction 
of  propagation.  Fortunately,  the  influence  of  lift  on  HSI  noise  appears  to  be  secondary22. 
Thus,  one  can  estimate  this  even  with  nonlifting  configurations. 

Recent  studies  using  the  Euler  and  Navier-Stokes  method  of  TURNS  have  demon¬ 
strated  the  feasibility  of  using  a  single  computational  fluid  dynamics  (CFD)  method  to 
calculate  simultaneously  the  aerodynamics  and  acoustics  of  a  helicopter  rotor  in  hover  and 
forward  flight23,  including  the  blade- vortex  interaction.  The  calculations  of  the  high-speed 
impulsive  noise  are  in  good  agreement  atleast  upto  a  distance  of  about  three  rotor  ra¬ 
dius.  This  numerical  procedure  is,  thus,  capable  of  calculating  the  aerodynamics  and  the 
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acoustics  (HSI)  out  to  several  rotor  diameters  all  in  one  single  solution  without  using  any 
modeling  either  for  wake  or  for  acoustic  propagation.  The  vortical  wake  and  its  influence  as 
well  as  the  acoustics  are  captured  as  a  part  of  the  overall  flowfield  solution.  Comparisons 
of  the  numerical  results  with  the  available  experimental  data  demonstrate  the  accuracy 
and  suitability  of  this  numerical  method.  The  next  few  sections  of  this  report  describe  the 
details  of  the  numerical  procedure  of  TURNS  code  and  example  solutions  are  given  for 
validation. 


2.  GOVERNING  EQUATIONS  AND  NUMERICAL  SCHEME 

The  governing  differential  equations  are  the  thin  layer  Navier- Stokes  equations.  These 
can  be  written  in  conservation-law  form  in  a  generalized  body- conforming  curvilinear  co¬ 
ordinate  system  as  follows24: 

9rQ  +  a(E  +  3,F  +  3(S=S.a(s  +  m  d) 

where  Re  is  the  Reynolds  number,  e  =  0  or  1  for  the  Euler  or  the  Navier-Stokes  equations, 
respectively,  and  r  =  t,  £  =  £(x,  y,  z,  t),  r\  =  r](x,  y,  r,  t)  and  C  =  ((*!  y>  ri  <)•  The  coordinate 
system  (x,y,  z,t)  is  attached  to  the  blade  (see  Fig.  1).  The  vector  of  conserved  quantities 
Q  and  the  inviscid  flux  vectors  E,  F ,  and  G  are  given  by 
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In  these  equations,  H  =  (e-f  p)  and  U,  V,  and  W  are  the  contravariant  velocity  components 
defined,  for  example,  as  U  —  4-  £xu  -i-  £yv  +  £zw.  The  Cartesian  velocity  components  are 

given  by  u,  v,  and  w  in  the  x,  y,  and  z  directions,  respectively.  Also,  the  density,  pressure, 
and  the  total  energy  per  unit  volume  axe  represented  by  p,  p,  and  e,  respectively.  While  the 
velocity  and  length  scales  are  nondimensionalized  by  the  characteristic  velocity  and  length 
scales,  given  by  the  ambient  sound  speed  a<x>  and  the  rotor  blade  chord  c,  the  pressure 
p,  density  p,  and  the  energy  e  axe  nondimensionalized  by  the  freest  ream  reference  values 
Pooh i  P°o,  and  pooa2oo,  respectively.  The  quantities  £z,  £*,  £t,  etc.  are  the  coordinate 

transformation  metrics  and  J  is  the  Jacobian  ofjthe  transformation.  For  the  thin  layer 
approximation  used  here,  the  viscous  flux  vector  S  is  given  by 
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where  Pr  is  the  Prandtl  number,  7  is  the  ratio  of  specific  heats,  and  a  is  the  speed  of  sound. 
For  the  noninertial  reference  frame  used  in  this  study,  source  terms  have  to  be  included 
in  Eq.  (1)  to  account  for  the  centrifugal  acceleration  of  the  rotating  blade25.  The  term 
represents  this  in  Eq.  (1)  and  is  given  for  a  rigid  rotor  rotating  in  the  x  -  y  plane  by 
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where  ft  is  the  angular  velocity  of  the  rotor.  The  fluid  pressure  p  is  related  to  the  conserved 
flow  quantities  through  the  nondimensional  equation  of  state  for  a  perfect  gas  given  by, 

p  =  (nf-l){e-~(u2  +  v2 +w2)}  (5) 

For  turbulent  viscous  flows,  the  nondimensional  viscosity  coefficient  p  in  S  is  computed  as 
the  sum  of  pi  +  pt  where  the  laminar  viscosity,  pi ,  is  estimated  using  Sutherland’s  law  and 
the  turbulent  viscosity,  pt,  is  evaluated  using  the  Baldwin- Lomax  algebraic  eddy  viscosity 
model26. 

A  finite-difference,  upwind  numerical  algorithm  is  developed  for  the  helicopter  rotor 
applications.  In  this,  the  evaluation  of  the  inviscid  fluxes  is  based  on  an  upwind-biased  flux- 
difference  scheme  originally  suggested  by  Roe27  and  later  extended  to  three-dimensional 
flows  by  Vatsa  et  al.20  The  chief  advantage  of  using  upwinding  is  that  it  eliminates  the  addi¬ 
tion  of  explicit  numerical  dissipation  and  has  been  demonstrated  to  produce  less  dissipative 
numerical  solutions20.  This  feature,  coupled  with  a  fine  grid  description  in  the  tip  region, 
increases  the  accuracy  of  the  wake  simulation.  The  van  Leer  monotone  upstream-centered 
scheme  for  the  conservative  laws  (MUSCL)  approach28  is  used  to  obtain  the  second-  or 
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third-order  accuracy  with  flux  limiters  so  as  to  be  total  variation  diminishing  (TVD).  Lower- 
Upper  -  Symmetric  Gauss-Seidel  (LU-SGS)  scheme,  suggested  by  Jameson  and  Yoon29-30, 
is  used  for  the  implicit  operator. 


The  space-discretized  form  of  the  differential  Eq.  (1)  at  node  (j,  k,  l)  is 


«  a  Ei±Ll£izi  Gt+j-Gt-j 

°rQ  A£  At)  AC 


(6) 


where  k,  and  /  correspond  to  the  C,  r/,  and  C  coordinate  directions,  respectively. 

The  application  of  Roe’s  upwinding  to  the  numerical  flux  of  the  inviscid  terms  results 
in  the  locally  one- dimensional  form  and  can  be  written,  e.g.,  in  the  f  direction,  as 


i[f(0„,(V{/J)/+J)  +  E(QL,(V(/J) li+t) 

-  -  ql)\  (7) 

where  A  is  the  Roe-averaged  Jacobian  matrix  and  Ql  and  Qr  are  the  left  and  right  state 
variables.  The  scheme  degenerates  to  the  first-order  accuracy  if  Ql  =  Q}  and  Qr  =  Qj+i, 
for  example,  at  the  grid  boundaries.  Higher-order  schemes  can  be  constructed  from  a  one- 
parameter  family  of  interpolations  for  the  primitive  variables  p,  u,  v,  w,  and  p.  For  example, 
the  left  and  right  state  variables  for  p  are, 

PL  =  {1  +  -J-K1  _  ")V  +  (1  +  (8“) 

PR  =  {1  -  ^p[(l  +  «)v  +  (1  -  *)A]}p,+l  (86) 

where  V  and  A  are  backward  and  forward  difference  operators,  and  k  is  a  parameter  that 
controls  the  construction  of  higher-order  differencing  schemes.  For  example,  k  =  j  is  used 
in  the  present  method  to  construct  the  third-order  scheme.  The  limiter  0  is  calculated  by 
using  Koren’s  differentiable  limiter31  as 

* _ 3Vp,&p,  +  ; _ 

‘  2(A Pi  -  Vp,f  +  3Vp;Ap,  +  <  W 


where  a  small  constant,  typically  e  =  10-6,  is  added  to  prevent  the  division  by  zero.  Similar 
formulae  axe  used  for  the  other  primitive  variables.  The  viscous  flux  terms  are  discretized 
using  standard  second-order  central-differencing24. 
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The  time  marching  integration  procedure  uses  the  LU-SGS  method.  The  details  of 
this  scheme  are  described  elsewhere21.  Briefly,  the  lower- upper  factored  symmetric  Gauss- 
Seidel  method  is  a  direct  modification  of  the  approximate  lower-diagonal-upper  (LDU) 
factorization  to  the  unfactored  implicit  matrix.  The  resulting  factorization  can  be  regarded 
as  the  symmetric  Gauss-Seidel  relaxation  method.  The  LDU  factorization  yields  better 
stability  than  the  simple  LU  factorization  since  the  diagonal  elements  always  have  the 
absolute  value  of  the  Jacobian  matrices. 

The  final  form  of  this  algorithm  can  be  written  for  a  first-order  time  accurate  scheme 
as 


LDUAQn  =  -A  tRHSn  (10) 

where 

L=I-  AtA~\J<kJ  +  AtV^A+ 

-AtB-\jXi  +  AtVvB+ 

—  AtC~\jtk,i  4-  A<V<C+  (11a) 

D  =(/  +  At(A+  -  A~ 

+  B+-B-+C+-C-)\j,k,i}~1  (116) 

U  —I  +  AtA+  Ij,*,/  ■+-  AtA^A 
+  Afi?+|Ji*i/  +  AtAvB~ 

^ ~  dR 

+  +  AtA^C  +  (He) 

where  At  is  the  time  step,  RHS  represents  the  discretized  steady  state  terms,  e.  gM  Eq.  (6), 
and  n  refers  to  the  current  time-level.  Also,  A+  =  1(A  +  cr^),  A~  =  1(A  —  <?(),  = 

\U\  +  ar $  -f-  e,  e  =  0.01  typically,  and  =  \J(X2  +  £y2  +  Also,  in  addition  to  the  source 
term  3?  added  to  the  right  side  of  the  Eq.  (6),  a  source  term  Jacobian,  dR/dQ ,  is  added  to 
either  L  or  U  operator  on  the  left  hand  side  of  Eq.  (10).  In  the  present  case  this  is  added 
to  the  U  operator,  as  shown  in  Eq.  (11c),  and  is  given  by 


■o 

0 

0 

0 

0- 

dR 

1 

0 

0 

Q 

0 

0 

0 

-Q 

0 

0 

0 

dQ 

J 

0 

0 

0 

0 

0 

.0 

0 

0 

0 

0. 

(12) 


As  a  result  of  the  simplified  form  of  the  Jacobian  terms,  e.g.,  A+,  all  the  diagonal  elements 
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of  L,  D,  and  U  reduce  to  scalar  elements.  Thus  this  method  requires  only  two  (one  forward 
and  one  backward)  sweeps  with  scalar  inversions  and  leads  to  less  factorization  error.  The 
source  term  adds  two  off-diagonal  elements,  resulting  in  slightly  more  computational  work. 

The  present  numerical  scheme  uses  a  modified  finite- volume  method  for  calculating  the 
metrics.  The  chief  advantage  of  standard  finite- volume  formulations  is  that  both  space  and 
time  metrics  can  be  formed  accurately25  and  that  the  free-stream  is  captured  accurately17. 
However,  to  be  compatible  with  the  present  finite-difference  numerical  scheme,  the  space 
metrics  are  evaluated  at  the  grid  nodes  instead  of  cell  interfaces.  But,  the  time  metrics  are 
evaluated  in  the  same  manner  as  in  a  finite- difference  scheme,  which  is  computationally  less 
expensive  than  a  rigorous  finite- volume  calculation.  As  a  result,  free-stream  subtraction  of 
the  time-metric  terms  is  then  required  to  restore  time  accuracy. 

In  the  calculation  procedure,  the  rotor  blade  is  started  from  rest  in  a  quiescent  fluid 
and  the  evolution  of  the  flowfield  is  monitored  as  the  blade  moves  in  azimuth.  To  take  ad¬ 
vantage  of  the  quasi-steady  nature  of  the  asymptotic  hovering  rotor  flowfield  in  blade-fixed 
coordinates,  a  locally-varying  time  step  is  used  in  the  integration  procedure  to  accelerate 
convergence,  as  suggested  in  Ref.  32. 

3.  GRIDS  AND  BOUNDARY  CONDITIONS 

Body-conforming,  single-block,  three-dimensional  computational  grids  were  con¬ 
structed  for  the  rectangular  rotor  blades33  by  stacking  and  bending  two-dimensional  C-H 
grids  which  were  generated  by  an  elliptic  solver34.  Because  of  the  cylindrical  nature  of 
the  flow  of  a  hovering  rotor,  a  C-H  cylindrical  grid  topology  was  chosen,  as  in  Ref.  17, 
with  the  wrap-around  C-direction  in  the  chordwise  direction  and  H-type  in  the  spanwise 
direction.  In  contrast  to  the  experimental  model  rotor  that  has  a  square  tip,  the  present 
grid  generator  approximates  this  as  a  bevel  tip  because  of  the  H- topology  in  the  spanwise 
direction  (see  Ref.  35). 

The  standard  viscous  grids  used  here  had  217  grid  points  in  the  wrap-around  (along  the 
chord)  direction  with  144  points  on  the  body,  71  points  in  the  spanwise  (radial)  direction 
with  55  points  on  the  blade  surface,  and  61  points  in  the  normal  direction.  The  grid  was 
clustered  near  the  leading  and  trailing  edges  and  near  the  tip  region  to  resolve  the  tip 
vortex.  It  was  also  clustered  in  the  normal  direction  to  resolve  the  viscous  flow  near  the 
blade  surface.  There  are  about  15  points  in  the  boundary  layer  with  a  spacing  of  the 
first  grid  point  from  the  surface  equal  to  5x  10-5  chord  (  y+  =  0(1)  ).  A  coarse  grid  was 
constructed  from  this  fine  viscous  grid  by  removing  every  other  point  in  all  three  directions. 
The  inboard  plane  near  the  axis  of  rotation  was  located  at  a  radial  station  equal  to  one 
chord.  The  grid  outer  boundaries  were  set  at  8  chords  in  all  directions.  The  same  grids 
were  used  for  the  Euler  calculations. 

Figure  1  shows  the  coarse  grid  that  was  used  in  the  computations.  Because  of  the 
symmetry  of  the  hovering  flow  and  the  periodic  boundary  condition  described  below,  the 
calculations  could  be  performed  for  only  one  blade.  Figure  la  shows  the  cylindrical  nature  of 
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the  grid  in  the  plane  of  the  rotor,  and  Fig.  lb  shows  an  isometric  view  of  the  grid  boundary 
for  a  single  blade.  For  clarity,  the  figure  shows  only  the  blade,  side  and  bottom  boundaries. 
Also  shown  is  the  coordinate  system,  where  x  is  in  the  chordwise  direction,  y  is  in  the  radial 
direction,  and  z  is  in  the  normal  direction.  The  blade  motion  is  counterclockwise. 

All  the  boundary  conditions  are  applied  explicitly.  At  the  wall  a  no-slip  boundary 
condition  is  used  for  the  viscous  calculations.  The  Euler  calculations  use  an  extrapolation 
of  the  contravariant  velocities  at  the  surface.  The  density  at  the  wall  is  determined  by 
a  zeroth-order  extrapolation.  The  pressure  along  the  body  surface  is  calculated  from  the 
normal  momentum  relation  (see,  for  example,  Ref.  24).  The  total  enery  is  then  determined 
from  the  equation  of  state.  To  ensure  continuity  across  the  wake  cut  and  also  outboard 
of  the  blade  tip,  where  the  grid  collapses  to  a  singular  plane  because  of  H-grid  topology, 
the  flow  quantities  are  determined  by  averaging  the  flow  variables  from  both  sides  of  the 
singular  plane. 

To  capture  the  information  in  the  wake  region  of  the  blade,  a  periodicity  condition  is 
used  in  the  blade  azimuthal  direction  that  swaps  the  flow  information,  after  interpolation, 
at  the  front  and  back  boundaries  of  the  cylindrical  grid  (see  Fig.  lb).  This  is  also  done  in 
an  explicit  manner.  The  flowfield  of  the  entire  rotor  is  then  assembled  by  combining  the 
flowfields  of  individual  blades  through  a  post-processing  procedure.  The  radial  inboard  and 
far-field  boundaries,  as  well  as  the  upper  boundary  of  the  cylindrical  mesh,  are  updated  by 
means  of  a  characteristic-type  boundary  condition  procedure,  although  the  Roe’s  upwinding 
used  in  the  numerical  procedure  would  otherwise  treat  the  boundaries  in  a  1-D  characteristic 
sense  anyway.  At  the  bottom  boundary,  the  scene  of  the  far-field  wake,  an  approximate 
condition  based  on  the  normal  velocity  is  used.  For  an  outflow  condition,  all  conserved  flow 
quantities  are  extrapolated  from  the  grid  interior  except  for  the  energy,  which  is  calculated 
by  prescribing  the  free-stream  pressure. 

Although  the  helicopter  rotor  operates  in  a  quiescent  fluid  atmospere,  unlike  in  a  fixed 
wing  airplane,  it  induces  significant  velocities  at  large  distances  from  the  rotor.  Therefore, 
specification  of  no  flow  at  the  inflow  boundaries  of  a  computational  box,  which  is  typically 
small  for  economy,  poses  a  difficulty  for  the  prediction  methods.  This  no-flow  condition  at 
the  farfield  boundaries  produces  a  “closed  box”  environment  for  the  rotor  where  the  flow 
recirculates  within  the  computational  “box”.  This  problem,  and  the  large  time  required 
for  the  initial  transients  to  decay,  were  recognized  by  Kramer  et  al.14,36,  who  used  an 
approximate  vortex-element  solution  to  specify  an  initial  condition  that  produced  a  flow 
through  the  farfield  boundaries  of  the  computational  box. 

A  simpler  and  more  economical  alternative  was  introduced  in  Ref.  37  using  the  con¬ 
cepts  of  a  three-dimensional  point-sink  and  simple  momentum  theory  as  a  guide.  This  is 
described  in  detail  later  in  the  text  in  Section  4.3.  With  this  approach,  the  application 
of  the  above  characteristic-type  boundary  condition  produced  a  non-zero  inflow  at  these 
boundaries.  At  the  outflow  boundary  located  at  the  farfield  boundary  below  the  rotor 
plane,  the  specification  of  the  flow  velocity  is  dictated  from  the  momentum  theory  con¬ 
cepts,  viz.,  the  flow  exit  through  a  circular  hole  whose  area  is  half  that  of  the  rotor  disk, 
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with  an  outflow  velocity  magnitude  twice  the  momentum  theory  average  value  through  the 
plane  of  the  rotor.  A  characteristic- type  numerical  outflow  boundary  condition  was  applied 
across  the  exit  plane  by  prescribing  this  outflow  mass  flux,  and  extrapolating  the  other  flow 
variables  from  within.  With  this,  the  flow  entered  and  exited  smoothly  the  computational 
box.  Most  of  the  results  presented  in  this  study  used  a  “box-type”  boundary  conditions  ap¬ 
proach  mentioned  earlier.  After  recognizing  the  importance  of  smooth  vortex  wake  descent 
below  the  rotor  disk,  one  calculation  was  repeated  with  the  new  boundary  conditions. 

4.  RESULTS  AND  DISCUSSION 

4.1  Model  Two-Bladed  Rotor  in  Hover 

The  test  cases  considered  in  this  study  correspond  to  the  experimental  model  hover 
test  conditions  of  Caradonna  and  Tung  (Ref.  33).  The  experimental  model  consists  of  a 
two-bladed  rigid  rotor  with  rectangular- planform  blades  with  no  twist  or  taper.  The  blades 
are  made  of  NACA  0012  airfoil  sections  with  an  aspect  ratio  of  6.  Three  experimental 
conditions  were  chosen  from  the  date:  1)  tip  Mach  number  Mup  =  0.44,  collective  pitch  9C 
=  8°,  and  the  Reynolds  number  based  on  the  blade  tip  speed  and  chord,  Re  =  1.92x10®; 
2)  Mtip  =  0.877,  9C  =  8°  and  Re  =  3.93x10®;  and  3)  M,ip  =  0.794,  9C  =  12°  and  Re  = 
3.55x10®. 

4.1.1  Fine  Grid  Results 

Surface  pressures  are  shown  in  Figs.  2-4  for  the  three  experimental  conditions  consid¬ 
ered.  These  calculations  were  done  on  a  fine  grid  consisting  of  nearly  one  million  points. 
Figure  2  shows  the  surface  pressures  for  the  conditions  of  Mup  =  0.44,  9C  =  8°,  and  Re 
=  1.92x10®.  In  this  figure,  the  present  calculations  are  compared  with  the  experimental 
data  of  Ref.  29  and  the  results  from  a  previous  Navier-Stokes  calculation  that  used  a  simple 
wake  model  (Ref.  9).  The  present  calculations  agree  well  with  the  experimental  data  for  all 
radial  stations.  There  are  some  improvements  in  the  results  at  y/R  =  0.50  and  0.96  over 
the  previous  results  from  Ref.  9.  It  should  be  pointed  out  that  the  calculations  of  Ref.  9 
used  a  0-0  grid  topology  with  nearly  700,000  grid  points  having  a  grid  clustering  similar 
to  the  present  grid. 

Figure  3  shows  a  comparison  of  surface  pressures  for  the  condition  of  M<ip  =  0.877,  9C 
=  8°  and  Re  =  3.93x10®.  At  this  transonic  flow  condition,  the  present  calculations  show 
excellent  agreement  with  the  experimented  data  for  all  radial  stations.  In  contrast  to  the 
calculations  of  Ref.  9,  the  present  results  show  shock  locations  and  shapes  that  are  well 
captured  due  to  the  TVD  upwinding  used  here.  The  inboard  regions  of  the  flow  are  also 
predicted  more  accurately;  this  indicates  that  the  present  computed  wake  is  superior  to  the 
approximate  wake  model  of  Ref.  9. 

Figure  4  shows  a  comparison  of  surface  pressures  for  the  condition  of  Mup  =  0.794,  9C 
=  12°  and  Re  =  3.55x10®.  Because  of  the  high  collective  pitch,  this  case  is  more  severe  in 
terms  of  the  shock  strength  and  shock-induced  boundary  layer  separation,  even  though  the 
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tip  speed  is  slightly  less  than  the  previous  case.  The  results  show  good  agreement  of  the 
calculations  with  the  experimental  data,  especially  near  the  tip. 

Figure  5  shows,  through  the  surface  particle  flow,  the  extent  of  shock-induced  boundary 
layer  separation  for  the  transonic  cases  discussed  above.  These  are  are  created  by  releasing 
fluid  particle  tracers  at  one  grid  point  above  the  surface  and  forcing  them  to  stay  in  that 
plane.  Figure  5a  shows  the  details  of  this  flow  for  the  case  of  M*,p  =  0.877  and  Bc  =  8°.  The 
separation  and  reattachment  locations  are  apparent  in  this  figure.  It  is  seen  that  this  flow 
condition  produces  a  mild  shock-induced  separation  in  the  outboard  part  of  the  blade.  In 
contrast,  the  shock-induced  separation  and  viscous-inviscid  interaction  are  much  stronger 
for  the  case  of  Mup  =  0.794  and  Bc  =  12°.  The  surface  particle  flow  pattern  for  this  more 
severe  case  is  shown  in  Fig.  5b.  As  seen,  the  extent  of  the  separation  is  much  larger  for 
this  flow  condition  than  for  the  case  of  Fig.  5a.  It  is  interesting,  however,  to  note  that  the 
separation  patterns  in  the  tip  region  are  approximately  the  same  for  these  cases. 

A  general  comparison  of  the  present  results  with  the  experimental  data  can  be  made  by 
examining  the  bound  circulation  distribution.  Figure  6  shows  such  a  plot  of  dimensionless 
circulation,  T/QR2,  as  function  of  r/R  for  6C  =  8°  case  and  tip  speeds  of  0.44  and  0.877, 
corresponding  to  the  data  presented  in  Figs.  2-3.  Here  r  is  the  radial  distance  from  the 
rotation  axis,  R  is  the  radius  of  the  rotor,  and  T  is  the  circulation,  computed  from  the  blade- 
element  lift.  Also  shown  are  the  integrated  data  from  the  experiments,  which  were  reported 
to  be  essentially  independent  of  tip  speed.  The  calculations  show  a  fair  agreement  with  the 
experiments,  except  in  the  inboard  part  of  the  blade.  This  suggests  that  only  the  near-field 
effects  of  the  tip  vortex  are  captured  as  well  as  desired.  There  are  two  possible  reasons  for 
the  poor  agreement  in  the  inboard  part  of  the  blade.  First,  the  vortex  wake  becomes  diffused 
in  the  far-field  grid,  so  its  induced  effect  is  significantly  diminished.  Second,  the  inboard 
plane  boundary  condition  may  be  indequate.  In  contrast  to  the  experimental  observation 
of  the  bound  cirulation  distribution,  the  present  calculations  show  some  dependency  on  the 
blade  tip  speed. 

In  the  tip  region  the  agreement  is  also  not  very  good.  This  may  be  due  to  the  bevel  tip 
that  is  used  in  the  computations  compared  to  the  square-tipped  blade  used  in  the  experi¬ 
ments.  Overall,  however,  the  surface  pressure  distributions  give  the  appearance  of  agreeing 
better  with  the  experiments  than  the  bound-circulation  distribution.  Relatively  minor 
discrepancies  in  the  pressure  distributions  near  the  leading  edge,  where  the  experimental 
transducer  locations  are  relatively  sparse,  seem  to  translate  into  significant  differences  in 
the  circulation  distribution. 

The  chief  advantage  of  the  Navier-Stokes  methods  is  to  predict  the  separated  flow  in 
the  tip  region  and  the  associated  detailed  structure  of  the  tip  vortex.  The  prediction  of 
the  overall  shed-wake  geometry  is  the  most  important  step  in  the  process  of  accurately 
modeling  the  complete  hover  flowfield.  The  ability  to  preserve  this  shed  wake  (including 
the  vortex  structure)  from  numerical  diffusion  is  a  more  complex  issue.  The  path  of  the 
tip  vortex  is  more  important  to  the  outboard  blade  loading,  whereas  the  ability  to  convect 
this  shed  wake  without  significant  numerical  diffusion  strongly  influences  the  inflow  in  the 
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inboard  parts  of  the  blade. 

Figure  7  shows  a  near-field  view  of  the  tip  vortex  particle  path  trajectory  for  the 
experimental  conditions  of  Mup  =  0.794  and  6C  =  12°  corresponding  to  Fig.  4.  These 
trajectories  Eire  generated  by  releasing  particle  tracers  in  the  vicinity  of  the  tip  of  the  blade 
on  both  surfaces  and  allowing  them  to  move  freely  in  time  and  space.  It  is  apparent  from 
this  figure  that  the  particles  released  right  on  the  tip  first  group  together  and  then  get 
braided  and  stay  in  the  vicinity  of  the  core.  As  observed  before  (Ref.  35),  the  process  of 
formation  of  the  tip  vortex  involves  braiding  of  fluid  particles  from  both  upper  and  lower 
surface  of  the  blade.  As  the  process  of  braiding  of  fluid  particles  from  upper  and  lower 
surfaces  continues,  the  tip  vortex  lifts  up  from  the  upper  surface  and  rolls  inboard  in  the 
downstream  wake. 

After  identifying  those  fluid  particles  released  in  the  vicinity  of  the  blade  tip  that  end 
up  in  the  vortex  core,  fewer  particles  were  released  on  the  tip  of  the  blade,  in  the  proximity 
of  the  quarter-chord  region,  to  trace  out  a  trajectory  of  the  tip  vortex  path.  The  free-wake 
trajectory  showed  the  correct  wake  contraction  and  descent  initially,  up  to  about  360  deg. 
of  vortex  age.  Subsequently,  the  wake  trajectory  continued  to  descend,  but  it  expanded  in 
size  and  eventually  ended  up  in  the  recirculating  flow.  This  problem,  due  to  the  “box-type” 
of  farfield  boundary  condition  used  here,  has  been  recently  corrected  by  applying  a  point- 
sink  and  simple  momentum  theory  concepts37,  as  noted  in  the  previous  Section.  With  this 
the  flow  smoothly  enters  the  computational  box  and  leaves  through  the  exit  boundary.  The 
calculated  vortex  trajectory  with  this  new  boundary  conditions  is  shown  in  Fig.  8.  The 
wake  contraction  and  descent  trajectories  are  shown  in  Fig.  8a.  Figure  8b  shows  a  very 
good  comparison  of  the  calculated  nearfield  trajectory  with  experiments33.  Comparison 
of  the  calculated  loads  with  the  two  methods  of  prescribing  farfield  boundary  conditions 
showed  only  a  5%  difference.  However,  the  wake  trajectory  had  the  right  characteristics  of 
contraction  and  descent  with  the  latter  approach  of  specifying  farfield  boundary  conditions. 

4.1.2  Fine  Grid  versus  Coarse  Grid  Results 

The  results  presented  in  the  preceeding  sections  were  calculated  on  a  fine  grid  of  nearly 
one  million  points.  The  initial  test  calculations  were  made  primarily  on  a  coarse  Navier- 
Stokes  grid  of  109x36x31  size.  This  grid  was  generated  by  removing  every  other  point 
from  the  fine  grid  in  all  three  directions.  The  outer  dimensions  of  the  grid  and  the  grid 
topology  are  thus  the  same  as  for  the  previous  fine  grid. 

Figure  9  shows  a  comparison  of  surface  pressure  distributions  for  the  fine  and  coarse 
grids  for  the  experiment eu  condition  of  Mup  =  0.877,  9C  =  8°  and  Re  =  3.93xl06.  A  general 
deterioration  of  the  predicted  surface  pressure  distributions  can  be  seen  for  the  coarse  grids. 
In  particular ,  the  shocks  are  not  as  sharp  as  for  the  fine  grid.  The  inboard  results,  not  shown 
here  for  y/R  <  0.5,  had  much  poorer  comparison.  The  tip  vortex  structure  is  also  very 
diffused  due  to  the  poor  grid  density  in  this  region.  However,  typical  converged  solutions 
with  this  coarse  grid  took  only  about  one  hour  of  CPU  (central  processor  unit)  time  on  the 
Cray-2  supercomputer,  down  from  about  15  hours  for  the  fine-grid  cases. 
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4.1.3  Euler  versus  Navier-Stokes  Results 


As  discussed  earlier,  there  have  been  several  attempts  to  capture  rotor  wakes  using 
Euler  methods  (Refs.  14-17).  The  vortex  formation  in  the  tip  region  of  a  wing  or  a  helicopter 
blade  is  a  result  of  complex  three-dimensional  separated  flow,  and  it  is  not  clear  how 
well  the  Euler  methods  are  able  to  mimic  viscous  effects  and  separation  to  produce  a 
vortex  structure.  Nevertheless,  these  Euler  methods  have  been  able  to  predict  the  pressure 
distributions  and  spanwise  blade  loading  reasonably  well  for  the  outer  part  of  the  blade. 
Against  this  background,  a  limited  comparison  of  surface  pressures  has  been  made  of  the 
Euler  results  with  experiments  and  Navier-Stokes  results.  The  Euler  calculations  were 
made  by  turning-off  viscous  terms  in  the  TURNS  code  and  using  the  same  fine  grid  of 
about  one  million  points  used  for  viscous  calculations.  It  should  be  noted  that  even  for  this 
fine  viscous  grid,  the  Euler  version  of  the  code  did  not  exhibit  any  stability  problems. 

A  typical  comparison  of  the  Euler  results  with  experiments33  and  the  Navier-Stokes 
results  is  presented  in  Fig.  10  for  the  test  conditions  of  Mtt>  =  0.877,  Bc  =  8°,  and  Re  = 
3.93xl06.  Because  it  neglects  viscous-inviscid  interaction,  the  Euler  method  overpredicts 
the  shock  wave  strength  and  position  (for  y/R  >  0.75).  Otherwise,  the  Euler  results 
are  in  good  agreement  with  the  Navier-Stokes  results,  which  show  mild  shock-induced 
separation  for  this  flow  condition  (see  Fig.  5a).  The  overall  agreement  of  surface  pressures 
certainly  does  not  reflect  the  details  of  the  flow  near  the  blade  surface,  especially  the 
separation  pattern  and  vortex  wake  details  as  predicted  by  the  Navier-Stokes  method.  A 
cursory  examination  of  surface  particle  tracers  showed  that  the  flow  is  completely  attached 
everywhere.  Although  it  is  not  clear  how  the  tip  vortex  is  formed  with  flow  separation  only 
at  the  edges  of  the  blade,  the  prediction  of  reasonably  accurate  surface  pressures  indicate 
that  the  induced  inflow  through  the  rotor  is  approximately  correct  and  that  the  tip  vortex 
strength  was  accurately  captured.  The  recent  Euler  calculations  of  Strawn38  indicate  that 
the  Euler  methods  capture  the  strength  and  peak  velocities  of  the  wing  tip  vortex  as  well  as 
any  Navier-Stokes  method  (e.  g.,  see  Ref.  35),  which  is  still  underpredicting  peak  velocities 
and  tight  braiding  of  the  vortex.  In  current  study,  the  details  of  the  wake  structure  for  the 
Euler  results  need  to  be  investigated  further. 

4.2  Four-Bladed  UH-60A  Blackhawk  and  BERP  Rotors  in  Hover 

Body  conforming  finite-difference  grids  are  used,  as  before,  for  the  long  slender  blades 
of  the  model  UH-60A  Blackhawk39  and  BERP40  rotors.  These  blades,  shown  in  planform 
in  Fig.  11,  have  aspect  ratios  of  15.51  and  12.46  respectively.  The  blades  are  highly  twisted 
and  have  twist  distributions  along  their  lengths  as  shown  in  Fig.  12.  Because  of  the  axial 
symmetry  of  the  hovering  rotor  flowfield,  a  C-H  grid  topology  is  chosen,  with  C-type  in  the 
chordwise  direction  and  H-type  in  the  spanwise  direction.  The  H-type  spanwise  topology 
creates  a  bevel  edge  for  the  square-tipped  UH-60A  blades.  To  adapt  to  the  cylindrical 
flowfield  of  the  rotor,  three-dimensional  grids  are  constructed  for  these  blades  by  stacking 
and  bending  two-dimensional  C-H  grids  so  that  they  lie  on  the  blade  along  the  chord. 
Typical  fine  grids  have  nearly  940,000  points  with  217  points  in  the  chordwise  direction 
and  71  and  61  points  in  the  spanwise  and  normal  directions,  respectively.  These  grids  have 


clustering  in  the  leading  and  trailing  edge  regions  as  well  as  in  the  tip  region.  The  grid 
lines  are  nearly  orthogonal  at  the  surface  with  a  spacing  of  0.00004  chord  in  the  normal 
direction  at  the  surface.  The  grid  boundary  is  located  at  two  rotor  diameters  away  at  the 
top  and  bottom  boundaries  and  one  radius  beyond  the  tip  in  the  spanwise  direction.  Coarse 
grids  have  been  constructed  by  eliminating  every  other  point  in  all  three-directions;  these 
typically  have  120,000  points.  Figure  13  shows  a  view  of  a  typical  grid  for  single  blade 
constructed  for  the  BERP  blade.  This  view  shows  the  grid  boundaries  and  the  details  of 
grid  clustering  near  the  grid  surface  and  in  the  tip  region. 

4.2.1  UH-60A  Blackhawk  Model  Rotor 

The  four-biaded  rotor  considered  here  is  the  model  Sikorsky  UH-60A  Black  Hawk 
rotor39.  The  model  blade  geometry  is  identical  to  the  full-scale  rotor  with  the  exception 
that  the  trim  tabs  have  been  omitted.  Therefore,  the  UH-60A  blades  used  here  have 
approximately  a  constant  chord  along  the  entire  span  as  shown  in  Fig.  1.  They  are  made 
up  of  SC1095  airfoils  except  for  0.48  <  r/R  <  0.87  where  SC1095R8  airfoils  are  used.  They 
have  a  swept-back  planform  in  the  tip  region.  A  distinguishing  feature  of  this  rotor  is  the 
hook-like  twist  distribution  in  the  tip  region  and  an  approximately  linear  twist  inboard  of 
r/R  ~  0.8,  as  shown  in  Fig.  12.  The  solidity  calculated  for  this  rotor  is  0.08393. 

At  the  time  these  calculations  were  performed,  much  of  the  experimental  information39 
for  the  UH-60A  model  rotor  was  restricted,  including  the  collective  pitch  settings  of  the 
model  rotor.  In  view  of  this,  a  representative  tip  speed  ( Mtip  —  0.628)  was  selected,  and  a 
thrust  condition  was  chosen  for  which  the  experimental  Figure  of  Merit  was  reported  to  be 
essentially  independent  of  thrust  coefficient,  i.e.,  near  FMma x-  A  collective  pitch  setting 
of  9C  =  9°  was  tried  out  to  match  with  the  experimental  data,  but  no  attempt  was  made 
to  model  the  elastic  twist  of  the  experimental  model  blades.  Figure  14  presents  the  calcu¬ 
lated  surface  pressure  distributions  at  several  radial  stations  along  with  the  experimental 
results  for  Ct/o  =  0.085,  which  approximately  matches  the  calculated  value  of  0.084.  The 
numerical  results  appear  to  be  in  general  agreement  with  the  data  at  all  radial  stations. 
The  normalized  sectional  thrust  distribution  along  the  blade  radius  is  also  compared  with 
the  experimental  data  in  Fig.  15.  The  overall  agreement  is  fairly  good.  The  calculated 
performance  parameters  for  this  rotor  are  presented  along  with  the  experimental  data  in 
Tables  1  and  2.  As  shown  in  Table  1,  the  calculated  Cq/<j  =  0.0068  agrees  well  with  the 
experimental  data. 

In  the  standard  definition,  the  total  power,  Cq,  is  the  sum  total  of  the  induced  and 
profile  powers.  It  is  also  the  sum  total  of  the  components  of  the  integrated  surface  pressures 
and  the  surface  skin  friction  Cq,  .  That  is,  Cq  =  Cq,  +  CQrro,  =  Cq9T  4-  Cq,  .  As  the 
integrated  surface  pressures  include  some  of  the  profile  power  also,  there  is,  of  course,  some 
degree  of  difficulty  in  isolating  the  actual  induced  power  and  profile  powers  according  to 
their  definitions.  Also,  the  grid  resolution  in  trailed  wake  precludes  accurate  determination 
of  induced  power  from  wake  integration,  but  the  induced  power  is  approximately  equal  to 
Cqvt  .  In  view  of  this,  the  second  definition  is  used  here  for  the  results  presented  in  these 
figures  and  in  Tables  1  and  2. 
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In  order  to  assess  the  importance  of  planform  effects,  flowfield  calculations  were  made 
on  an  equivalent  rectangular  blade  of  identical  airfoils  and  twist  distribution  as  the  UH- 
60A  blade  (shown  as  Twistl  in  Fig.  12)  for  the  above  flow  conditions  of  Mtip  =  0.628, 
and  8C  =  9°  and  Re  =  2.75  x  106.  The  results  are  compared  with  the  UH-60A  blade  in 
Fig.  15  in  a  plot  of  normalized  sectional  thrust  coefficient  (C*)  variation  along  the  blade 
radius  (r/R).  Since  the  equivalent  rectangular  blade  has  the  identical  twist  distribution  as 
the  UH-60A  blade  along  its  radius  and  has  an  unswept  planform  in  the  tip  region,  the 
differences  seen  with  the  UH-60A  blade  in  Fig.  15  axe  only  in  this  region.  The  sectional 
thrust  for  the  equivalent  blade  has  a  smaller  peak  than  the  UH-60A  blade  and  this  peak 
occurs  towards  the  tip  of  the  blade  as  seen  in  Fig.  15.  Although  both  of  these  blade  tips 
are  modeled  by  bevelled  edges,  because  of  the  H-grid  topology  used  in  the  radial  direction, 
the  UH-60A  blade  produces  a  relatively  milder  separation  in  this  region  than  the  equivalent 
rectangular  blade  with  am  unswept  tip.  Examination  of  the  flow  particle  tracers  released 
in  the  tip  region  for  these  blades  showed  that  the  tip  vortex  location  for  the  UH-60A  blade 
is  slightly  inboard  of  that  seen  for  the  equivalent  rectangular  blade.  This  observation  is  in 
general  agreement  with  some  previous  studies  of  swept  and  unswept  planforms  with  bevel- 
tips33.  In  terms  of  overall  performance,  the  equivalent  rectangular  blade  produces  about 
7%  less  thrust  than  the  UH-60A  blade,  and  the  calculated  Figure  of  Merit  is  about  5% 
lower.  Tables  1  and  2  detail  the  performance  parameters  for  these  rotors. 

Figure  16  shows  the  sectional  thrust  distribution  along  the  radius  of  UH-60A  blade  for 
two  thrust  conditions  corresponding  to  6C  =  9°  and  13°  at  Mtip  =  0.628  and  Re  =  2.75  x  106. 
The  higher  collective  pitch  produces  higher  thrust  all  along  the  blade  radius  and  in  all  nearly 
67%  more  thrust  than  the  8C  =  9°  case,  but  the  Figure  of  Merit  is  almost  unchanged.  The 
performance  parameters  for  this  thrust  condition  are  also  summarized  in  Tables  1  and  2. 
This  high  thrust  condition  also  produces  a  much  stronger  tip  vortex  which  is  located  more 
inboard  than  the  8C  =  9°  case.  The  flow  topology  is  very  complicated  for  the  upper  surface 
of  the  blade  at  this  high  thrust  condition  with  shock-induced  boundary  layer  separation 
occuring  over  a  large  part  of  the  blade,  as  shown  by  the  surface  particle  tracers  (skin  friction 
lines)  in  Fig.  17.  In  contrast,  as  seen  in  the  close-up  view  of  the  tip  region  in  Fig.  18,  the 
low  thrust  case  (8C  =9°)  does  not  produce  shock-induced  flow  separation.  Both  thrust 
conditions  produce  flow  separation  on  the  inboard  parts  of  the  blade  as  seen  in  Fig.  17, 
but  the  spanwise  extent  is,  of  course,  greater  at  8C  =  13°.  An  examination  of  the  Mach 
contours,  on  the  upper  surface  of  the  blade  at  a  plane  outside  the  boundary  layer  as  seen 
in  the  views  of  Fig.  19,  for  the  8C  —  13°  case  shows  the  extent  and  severity  of  the  transonic 
flow  with  shock  wave  and  shock  induced  separated  flow  around  the  swept  part  of  the  blade. 

Two  representative  performance  curves  for  the  UH-60A  rotor  as  calculated  here  for 
Mtip  =  0.628  are  shown  in  Fig.  20  and  Fig.  21  along  with  the  experimental  data6.  The 
details  of  the  calculated  results  are  as  shown  in  Tables  1  and  2.  Figure  20  shows  a  plot 
of  power  vs  thrust  for  this  rotor.  Also  shown  in  this  are  the  induced  power  calculated 
using  an  empirical  relation41  and  the  ideal  power.  The  present  calculations  for  thrust  and 
total  power  at  8C  =  9°  show  good  agreement  with  experiments.  Also,  at  this  subcritical 
flow  condition,  the  component  of  the  power  calculated  from  the  pressure  integration,  CQrr , 
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is  in  good  agreement  with  the  empirically-estimated  induced  power41.  The  high  thrust 
condition  (8C  =  13°)  for  this  rotor  is  apparently  beyond  the  model  test  range.  This  high 
thrust  condition  produces  severe  shock-induced  flow  separation  over  a  large  part  of  the 
blade  and  a  more  severe  flow  separation  pattern  on  the  inboard  part  of  the  blade.  As  seen 
in  Table  2,  the  skin  friction  contribution  to  the  power  has  slightly  decreased  due  to  the 
reversed  flow  in  the  separated  regions. 

Figure  21  shows  another  rotor  performance  plot  of  Figure  of  Merit.  As  before,  at  the 
moderate  thrust  condition  corresponding  to  9C  —  9°,  the  calculated  Figure  of  Merit  is  in 
good  agreement  with  experiments,  including  the  empirical  estimate  of  Figure  of  Merit  that 
neglects  viscous  contribution  ( FMi ).  It  is  interesting  to  note  that  the  calculated  Figure 
of  Merit  including  viscous  effects  is  almost  the  same  at  9C  =  13°  as  at  6C  =  9°,  indicating 
that  this  particular  measure  of  total  hover  performance  seems  to  stay  relatively  uniform  for 
this  rotor,  at  least  over  this  range  of  0C.  On  the  other  hand,  the  Figure  of  Merit  calculated 
using  only  Cg,,  from  Table  2  shows  a  significantly  lower  value  at  this  high  thrust  condition 
than  at  Bc  =9°.  In  other  words,  the  calculations  confirm  that  the  relative  contribution  of 
skin  friction  drag  to  Figure  of  Merit  is  greater  at  lower  thrust  conditions. 

In  addition  to  the  planform  change  study  for  this  blade,  a  fictitious  UH-60A  blade 
having  the  same  airfoils  and  planform  but  with  a  linear  twist  distribution  in  the  tip  region 
(Twist2  in  Fig.  12)  was  considered  to  examine  the  importance  of  twist  in  this  region  on 
airloads.  The  sectional  thrust  distribution  for  this  blade  is  shown  in  Fig.  15  along  with 
results  for  the  standard  UH-60A  blade  and  the  equivalent  rectangular  UH-60A  blade  at 
the  same  test  conditions  as  of  Fig.  14.  As  seen  in  this  figure,  this  blade  generates  almost 
identical  thrust  in  the  inboard  region.  There  appears  to  be  a  small  difference  in  the  tip 
region,  resulting  in  slightly  less  overall  thrust  coeficient  than  the  UH-60A  blade  and  with 
a  7%  lower  Figure  of  Merit,  as  shown  in  Table  1.  As  before,  there  is  no  flow  separation 
on  the  blade  except  in  the  very  tip  region  at  this  flow  condition.  An  inspection  of  particle 
tracers  released  in  the  tip  region  from  both  sides  of  the  blade  indicates  the  location  of  tip 
vortex  for  this  configuration  to  be  slightly  inboard  of  the  UH-60A  blade  and  has  less  tight 
braiding  than  the  UH-60A  blade  tip  vortex  at  a  comparable  blade  location,  although  this 
blade  also  has  a  swept-tip  planform. 

4-2.2  BERP  Rotor 

The  BERP  rotor  has  an  unorthodox  planform  shape  with  high-performance  RAE  air- 
foils.The  spanwise  twist  distribution  used  in  the  present  investigation  and  in  Ref.  40  is 
shown  in  Fig.  12;  it  corresponds  to  a  version  of  the  BERP  blade  designed  for  good  hover 
performance,  rather  than  the  version  used  on  the  Westland  G-Lynx  helicopter42,  and  the 
total  twist  is  approximately  the  same  as  the  UH-60A.  However,  the  innovative  distribution 
of  advanced  RAE  airfoils  is  used  for  the  present  BERP  blade.  As  noted  in  Ref.  40,  a  dis¬ 
tinguishing  feature  of  the  BERP  blades  is  the  enlarged-chord  paddle-shaped  planform,  (see 
Fig.  11),  with  a  leading-edge  notch  at  approximately  80%  span  and  a  swept  delta- wing-like 
planform  in  the  tip  region.  The  combination  of  unusual  planform,  high  performance  air- 
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foils,  and  a  balanced,  spanwise  distribution  of  airfoils  has  enabled  the  low-twist  BERP  rotor 
to  achieve  high  performance  in  high-speed  forward  flight.  In  order  to  compare  planform 
effects  of  the  high-twist  BERP  blade  with  the  UH-60A  blade,  two  calculations  are  made 
for  collective  pitch  of  6C  =  9°  and  13°  at  a  tip  speed  of  Mt,p  =  0.628  and  Re  =  2.95  x  106. 
The  calculated  rotor  solidity  is  0.10913. 

Figure  22  shows  sample  surface  pressures  for  the  BERP  rotor  at  different  radial  loca¬ 
tions  for  the  flow  condition  corresponding  to  9C  =  13°.  Currently  no  experimental  data 
are  available  for  these  flow  conditions.  However,  these  results  and  the  surface  pressures 
for  $c  =  9°  appear  to  be  in  good  qualitative  agreement  with  the  nonrotating  blade  results 
of  Ref.  40.  The  spanwise  circulation  distribution  is  shown  for  this  blade  for  6C  =  9®  in 
Fig.  23,  along  with  UH-60A  results  at  identical  collective  pitch  and  flow  conditions.  Since 
the  BERP  blade  has  a  variable  chord  along  its  radius  in  contrast  to  the  UH-60A  blade, 
a  comparison  of  bound  circulation  distribution  appears  to  be  more  meaningful  than  the 
sectional  thrust  distribution.  As  seen  in  the  comparison  of  Fig.  23,  the  combination  of  high 
performace  airfoils,  the  leading-edge  notch  and  increased  area  of  the  paddle  part  produces 
an  increased  thrust  on  most  of  the  blade  except  for  r/R  <  0.4  where  it  is  nearly  the  same 
as  the  UH-60A  blade.  The  increased  circulation  seen  at  the  inboard  station  (Ri)  may  be 
due  to  the  deficiency  in  the  inboard  boundary  condition. 

Tables  1  and  2  summarize  important  hover  performance  results  for  this  rotor.  Fig¬ 
ures  24  and  25  compare  the  sectional  thrust  and  bound  circulation  distributions  for  the 
two  thrust  conditions.  The  high  thrust  condition  produces  about  64%  more  thrust  at  the 
expense  of  101%  more  power  compared  to  the  moderate  thrust  condition  of  9C  =  9°.  The 
sectional  thrust  distributions,  shown  in  Fig.  24,  has  the  same  qualitative  behavior  as  that  of 
UH-60A  blade.  The  bound  circulation  distribution  of  Fig.  25  appears  to  resemble  qualita¬ 
tively  the  sectional  lift  distribution  of  nonrotating  BERP  blade40.  Of  particular  significance 
is  the  distribution  in  the  vicinity  of  the  notch.  As  observed  for  a  nonrotating  blade  from 
both  the  wind  tunnel  tests  and  calculations40,  the  present  calculations  also  show  a  decrease 
and  then  an  increase  of  sectional  thrust  distribution  in  this  region  for  the  dc  =  13°  case.  In 
contrast,  the  thrust  distribution  for  6C  —  9°  case  does  not  show  this  behavior. 

The  surface  particle  tracer  pictures  of  the  BERP  blade  are  shown  in  Fig.  26  for  6C  =  9° 
and  13°.  For  both  thrust  conditions,  the  flow  stays  attached  on  the  paddle  part  of  the  blade. 
As  seen  in  Fig.  26(a),  the  flow  appears  to  be  completely  attached  at  the  moderate  thrust 
condition.  The  high  thrust  condition  corresponding  to  9C  =  13°,  however,  produces  a 
complicated  vortical  flow  structure  in  the  notch  region  as  shown  in  Fig.  26(b)  as  well  as 
in  Fig.  27.  Also,  the  strong  supercritical  flow  at  this  high  thrust  condition  produces  a 
relatively  small  extent  of  shock-induced  flow  separation,  confined  only  to  the  inboard  of  the 
notch  region  as  seen  in  Fig.  15(b)  and  in  Figs.  27(a)  and  27(b).  In  contrast,  as  shown  in 
Fig.  18(b),  the  UH-60A  blade  produces  a  severe  shock  induced  flow  separation  over  a  large 
part  of  the  blade  for  this  flow  condition.  The  extent  and  severity  of  the  supercritical  flow 
for  the  upper  surface  of  the  BERP  blade  in  the  notch  region  is  highlighted  in  Fig.  27(c)  by 
the  pressure  contour  plot. 
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The  nearly  total  absence  of  the  flow  separation  on  the  paddle  part  of  the  BERP  blade  at 
the  transonic  flow  condition,  except  in  the  very  tip  region  where  the  tip  vortex  develops,  has 
a  favorable  influence  on  the  performance  of  the  blade.  Since  these  blades  have  high  twists, 
both  thrust  conditions  produce  severe  flow  separation  inboard.  The  extent  of  separation 
for  the  9C  =  13°  case,  shown  in  Fig.  27(a),  is  much  more  massive  than  the  9C  =  9°  case,  as 
for  the  UH-60A  blade  (see  Fig.  17). 

A  detailed  comparison  of  surface  particle  traces  for  the  UH-60A  (Fig.  18)  and  BERP 
(Fig.  26)  rotors  in  the  outer  part  of  the  blades  towards  the  tip  reveals  several  noteworthy 
features,  resulting  from  the  different  planforms  and  twist  distributions  in  this  region.  The 
UH-60A  blade  has  a  hook-like  twist  distribution  in  the  swept  part  near  the  tip  while  the 
BERP  blade  has  constant  twist  with  a  large  chord  in  the  paddle  part  and  a  delta-wing  like 
tapered  tip.  For  the  low  thrust  case,  corresponding  to  9C  =  9°,  both  blades  show  completely 
attached  flow  except  for  the  extreme  tip  region,  and  benign  surface  flow.  In  contrast,  at 
the  high  thrust  condition  (0C  =  13°),  they  have  totally  different  surface  flow  patterns.  The 
UH-60A  blade  shows  strong  shock-induced  separation  and  reattachment  near  the  leading 
edge  area  on  either  side  of  the  beginning  of  the  sweep-back,  whereas  the  BERP  blade  shows 
completely  attached  flow  on  the  paddle-part  except  in  the  extreme  tip  region.  The  notch 
region  has  a  complicated  flow  separation  and  vortical  flow  structure. 

Figure  28  shows  the  details  of  the  tip-vortex  formation  and  roll-up  process  for  the 
BERP  blade.  This  flow  visualization  picture  is  generated  by  releasing  unrestricted  flow 
particle  tracers  in  the  tip  region  on  both  sides  of  the  blade  at  different  chordwise  and 
spanwise  locations  and  at  different  heights  from  the  surface  of  the  blade.  The  particles  so 
released  bunch-up  and  get  braided  over  each  other  to  form  a  distinct  structure  of  concen¬ 
trated  vorticity.  Such  a  structure  represents  the  formation  process  of  a  concentrated  tip 
vortex  as  seen  in  Fig.  28.  While  this  tip  vortex  is  still  forming,  it  starts  to  roll-up  and  in  the 
process  stays  distinctly  above  the  wake  vortex  sheet.  Once  it  leaves  the  surface,  it  convects 
with  the  wake.  As  this  vortex  wake  convects,  it  contracts  in  size  and  descends  with  the 
flow.  With  the  present  nonadaptive,  structured  grids,  the  computed  vortex  structure  is 
smeared,  meaning  the  core  of  such  a  vortex  is  diffuse.  But  a  line  integration  of  the  velocity 
vector  along  a  closed  path  around  such  a  line  vortex  reveals  that  the  circulation  is  correctly 
captured  in  strength.  This  estimation  has  been  verified  at  several  locations  along  the  vortex 
wake  trajectory,  viz.,  in  the  near  wake  and  near  the  periodic  boundaries.  Because  of  the 
vortex  diffusion,  the  farther  the  integration  station  from  the  blade  the  larger  the  size  of 
the  closed  path  is  needed.  In  comparing  the  vortical  flow  structures  of  UH-60A  and  BERP 
blades  corresponding  to  9C  =  13°,  it  was  found  that  both  have  distinct  but  diffuse  tip  vortex 
structures.  The  notch  area  of  the  BERP  blade  produces  a  complex  vortical  structure,  which 
is  barely  discemable  in  Fig.  28.  Similary,  the  kink  of  the  UH-60A  blade  also  produces  a 
very  mild  vortical  structure  for  this  flow  condition. 

4.3  Effect  of  Farfield  Boundary  Conditions 

A  difficult  aspect  of  numerical  prediction  of  hovering  rotor  flows  is  the  prescription  of 
farfield  boundary  conditions.  In  hover,  if  the  specification  of  farfield  boundary  conditions 
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is  based  on  the  quiescent  flow  outside  of  the  computational  box,  then  the  flow  into  and  out 
of  the  box  will  be  zero.  This  creates  a  closed  box  environment  for  the  rotor,  with  the  flow 
inside  the  box  recirculating  when  the  rotor  is  spinning.  In  principle,  this  is  precisely  the 
environment  in  many  hover  test  chambers.  However,  a  large  computational  box  of  this  size 
for  the  Navier-Stokes  calculations  is  not  cost  effective  and  exceeds  the  memory  limits  of  the 
present  day  supercomputers.  Alternatively,  with  a  smaller  computational  box  around  the 
rotor  one  should  specify  a  non-zero  flow  farfield  boundary  conditions  that  will  allow  the 
flow  to  enter  and  exit  the  box  without  violating  the  conservation  laws.  Currently,  no  precise 
methods  are  available  to  specify  such  boundary  conditions,  although  some  schemes  have 
been  suggested43.  The  present  calculations  take  a  slighty  different  approach  to  specifying 
farfield  conditions.  A  three-dimensional  point  sink,  whose  strength  is  determined  from  the 
thrust  of  the  rotor  by  satisfying  the  mass  flow  conservation,  is  located  at  the  axis  of  rotation 
in  the  plane  of  the  rotor.  The  point  sink  pulls  the  flow  from  outside  into  the  computational 
box  with  a  velocity  W,n  given  by 


where  cP  —  x2  +  y2  +  z2.  The  velocity  given  by  Eq.  (13)  represents  the  magnitude  of  the  to¬ 
tal  incoming  flow  volocity.  An  appropriate  component  of  this  flow  enters  the  computational 
box  from  each  boundary  except  through  an  exit  plane  on  the  lower  boundary.  Assuming 
the  far-field  exit  velocity  at  this  exit  plane  to  be  uniform,  its  magnitude  can  be  determined 
from  1-D  momentum  theory  concepts  by  relating  the  outflow  momentum  to  the  rotor  thrust 

by 


A  characteristic- type  outflow  boundary  condition  is  specified  by  prescribing  this  w- velocity, 
We ,  across  the  exit  plane  whose  area  is  half  of  the  rotor  disk.  The  other  four  quantities  are 
extrapolated  from  within.  Prescription  of  this  new  boundary  condition  procedure  enables 
the  flow  to  enter  and  exit  smoothly  the  computational  box.  A  schematic  of  this  new  farfield 
boundary  of  the  computational  box  is  shown  in  Fig.  29. 

Using  the  above  approach  of  specifying  the  farfield  boundary  conditions,  calculations 
were  repeated  for  the  UH-60A  blade  at  flow  conditions  of  Mup  =  0.628,  9C  =  9°  and 
Re  =  2.75  x  106.  Figure  30  shows  the  sectional  thrust  distribution  calculated  with  the 
new  farfield  boundary  condition  compared  with  old  boundary  conditions  results  along  with 
experiments.  Although  there  is  general  agreement  of  the  overall  results  in  this  plot,  there 
are  small  differences  near  the  kink  and  in  the  tip  region  that  appear  to  be  due  to  differences 
in  the  wake  trajectory.  The  performance  parameters  for  this  new  condition  are  shown  in 
Table  1.  The  thrust  and  power  values  Eire  about  3%  more  than  those  obtained  with  the  old 
boundary  conditions. 


The  important  test  of  this  new  boundary  condition  is  the  behavior  of  the  captured 
wake  trajectory.  Fluid  particle  tracers  that  reside  within  the  core  of  the  tip  vortex,  such 
as  shown  in  Fig.  28,  were  identified  and  monitored  as  the  wake  evolved.  Figure  31  shows 
the  views  of  vortical  wake  contraction  and  descent  calculated  with  this  new  procedure. 
The  nearfield  trajectory  of  this  wake,  for  both  contraction  and  descent,  appears  to  be  in 
reasonable  agreement  with  the  experiments,  as  shown  in  Fig.  31c.  The  previous  calculations, 
for  the  same  computational  box  with  old  boundary  conditions,  produced  a  wake  trajectory 
that  was  accurate  only  up  to  about  180“  of  vortex  age  (azimuthal  travel).  Subsequently, 
the  wake  trajectory  continued  to  descend  but  expanded  in  size,  from  its  previous  size  at 
180°  azimuth,  and  eventually  got  caught  in  the  recirculating  flow.  With  the  new  farfield 
boundary  conditions,  the  entire  wake  exits  smoothly  at  the  bottom  boundary,  as  seen  in 
Fig.  31. 

4.4  Model  Two-Bladed  Rotor  in  Forward  Flight 

In  this  section,  instantaneous  surface  pressures  of  a  two-bladed  rotor  in  forward  flight, 
calculated  using  the  TURNS  code,  are  presented  and  compared  with  experimental  data 
for  rotor-alone  configuration44,45.  Several  test  conditions  from  among  many  model  rotor 
experiments  are  selected  for  computations.  The  hovering  rotor  calculations  described  in 
the  previous  sections  were  done  in  the  blade-fixed  coordinates  system.  In  contrast,  the 
forward  flight  calculations  are  done  in  the  inertial  coordinate  system.  The  TURNS  code 
switches  from  one-system  to  the  other  by  choosing  the  proper  input  for  the  grid  motion.  To 
validate  the  TURNS  code  for  forward  flight  motion,  calculations  are  done  for  a  nonlifting 
rotor  in  forward  flight  at  both  subcritical  and  supercritical  flow  conditions.  This  solution 
is  then  used  as  the  baseline  solution  for  calculating  a  lifting,  three-dimensional,  unsteady, 
blade-vortex  interaction  (BVI)  flowfield. 

4.4.1  Rotor  in  Forward  Flight 

The  two-bladed  rigid  rotor  considered  here  for  flowfield  calculation  is  a  nonlifting  con¬ 
figuration  that  was  tested  in  a  wind  tunnel  by  Caradonnaet  al.  44 .  Several  test  calculations 
have  been  performed46,47  corresponding  to  the  model  rotor  test  conditions44,45.  Typical 
results  of  instantaneous  surface  pressures  for  rotor- alone  case  are  presented  in  Figs.  32  and 
33  for  the  conditions  of  Mtip  =  0.6  and  0.8,  respectively,  and  an  advance  ratio  n  =  0.2. 
The  results  for  subcritical  case  of  Fig.  32  shows  that  the  rotor  flowfield  behaves  like  a 
two-dimensional  flow  and  the  surface  pressures  at  different  azimuthal  location  of  the  rotor 
all  look  alike.  The  unsteady  time-lag  effects  and  the  three-dimensional  effects  are  nearly 
absent.  In  contrast,  the  surface  pressures  shown  for  various  azimuthal  locations  of  rotor  in 
Fig.  33,  for  the  supercritical  flow  condition,  appears  to  behave  differently.  The  presence  of 
transonic  shock  wave  on  the  advancing  cycle  produces  strong  time-lag  effects  and  accentu¬ 
ates  three-dimensionality  of  the  flow  (see  also  Ref.  46).  The  instantaneous  surface  pressures 
for  subcritical  and  superctical  cases  are  in  good  agreement  with  experimental  data.  The 
numerical  method  captures  correctly  the  evolution  of  the  unsteady  flow  in  forward  flight, 
atleast  for  this  nonlifing  configuration. 
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4.4.2  Three-Dimensional  Parallel  Blade- Vortex  Interaction 


Using  this  baseline  solution,  a  vortex  fitting  technique48,49  is  used  to  introduce  the 
convecting  vortex  flowfield  in  to  the  rotor  flowfield  and  the  evolution  of  the  interactional 
flowfield  is  then  calculated.  The  interacting  line  vortex  lies  along  y-axis  and  is  located  at 
x0  =  0.0  and  z0  =  -0.4  where  x,  y,  and  z  are  the  physical  space  coordinates.  These  are 
respectively  along  the  chord  of  the  blade,  along  the  radius  of  the  blade  and  normal  to  the 
blade  surface.  Both  parallel  and  oblique  blade- vortex  interactions  have  been  calculated  for 
several  tip  Mach  numbers.  In  contrast  to  the  subcritical  flow  conditions,  the  presence  of 
shock  waves  under  supercritical  flow  conditions  appears  to  enhance  the  flow  unsteadiness 
and  three-dimensionality.  A  typical  result  of  instantaneous  surface  pressure  distribution  for 
a  parallel  blade- vortex  interaction,  calculated  on  a  coarse  grid,  is  shown  in  Fig.  34  along  with 
comparison  to  experimental  data.  For  the  azimuthal  position  of  the  blade  of  rp  =  174.5°, 
the  blade-vortex  interaction  appears  to  have  peaked.  As  seen,  the  instantaneous  surface 
pressures  are  in  fair  agreement  with  experiments.  Although  not  shown  here,  the  oblique 
BVI  occurs  over  a  larger  azimuthal  travel  of  the  blade.  The  interaction  begins  in  the  first 
quadrant  and  completes  only  towards  the  end  of  second  quadrant. 

4.5  High-Speed  Impulsive  Noise  in  Hover  and  Forward  Flight 

The  capability  of  TURNS  code  in  simulataneously  calculating  the  aerodynamics  and 
aeroacoustics  in  one  single  calculation  is  demonstrated  here.  The  acoustics  calculations 
were  performed  by  Baeder23.  For  calculating  the  HSI  noise,  the  rotor  blade  considered 
is  from  the  experiments  of  Ref.  50.  The  blade  is  a  l/7th  scale  model  of  a  UH-1H  main 
rotor  with  untwisted  rectangular  blades  of  NACA  0012  airfoil  section.  This  rotor  has  an 
aspect  ratio  of  13.7.  The  computations  were  done  on  a  49x37x55  grid  (half  plane)  with  grid 
boundary  located  at  48  chords  away  from  the  center  of  rotation.  The  first  case  examined 
is  for  Mtip  =  0.90  and  Re  =  1.6x10®  for  which  the  flow  is  just  delocalized.  The  pressure 
time  histories  at  3.09  rotor  radii  are  shown  in  Fig.  35a.  The  predicted  amplitude  of  the 
negative  peak  pressure,  as  well  as  the  wave  shape,  agree  well  with  the  experiments.  Similar 
agreement  occurs  at  other  locations  also.  For  this  weak  shock  case,  the  Euler  and  Navier- 
Stokes  results  are  nearly  identical.  In  contrast  to  the  present  CFD  method,  studies  using 
the  acoustic  analogy  method22  or  the  nonlinear  Kirchoff  formulation51,  with  aerodynamic 
input  from  either  small  disturbance  code  or  full  potential  code,  underpredict  the  negative 
peak  pressure. 

If  Mtip  is  increased  to  0.95,  the  shock  on  the  surface  of  the  blade  becomes  very  strong. 
Examination  of  the  pressure  contours,  in  a  plane  normal  to  the  blade  surface  at  98%  of  the 
span,  indicates  that  the  Navier-Stokes  results  show  the  shock  location  slightly  upstream 
and  more  smeared  at  its  root  them  the  Euler  results.  However,  away  from  the  blade  surface 
the  shock  location  and  strength  appears  to  be  nearly  identical  with  that  of  Euler  results. 
The  resulting  pressure  time  histories  at  3.09  rotor  radii  are  shown  in  Fig.  35b.  Again 
the  agreement  with  experiment  is  excellent.  The  reason  for  the  overall  good  agreement 
of  the  Euler  results  with  the  Navier-Stokes  results  is  because  HSI  noise  is  caused  by  the 
(summation  of  the)  large  gradients  across  the  whole  flowfield,  and  not  just  the  portion  of 
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the  shock  that  intersects  the  blade  surface. 


For  calculating  the  HSI  find  its  propagation  in  forward  flight,  the  OLS  rotor  blade  of 
an  AH-1  helicopter  is  selected  due  to  its  simple  geometry  of  a  rectangular  planform  with 
constant  thickness.  The  maximum  thickness  of  the  airfoil  is  9.71%  chord  and  the  rotor 
blade  has  an  aspect  ratio  of  9.22.  Acoustic  data  are  available  for  full-scale  and  model  OLS 
blades52.  For  the  calcuations  reported  here  the  twist  of  the  blade  is  neglected  and  the  tip 
path  plane  angle  and  collective  pitch  are  set  to  zero.  The  mesh  consists  of  73x55x55  points 
with  55x19  points  on  the  blade,  covering  the  upper  half  plane  of  the  whole  flowfield. 

The  test  case  examined  is  for  Mup  =  0.665,  an  advance  ratio  fi  =  0.348,  and 
Re  =  2.17xl06.  The  flowfield  is  delocalized22  for  a  short  period  of  time.  The  time  his¬ 
tories  of  the  pressure  pulse  at  1.8  rotor  diameters  directly  in  front  of  the  blade  are  shown 
in  Fig.  36.  Non-linear  effects  manifest  themselves  in  the  steeper  compression  wave  and 
in  the  magnitude  of  the  pulse;  calculations  using  linear  theory  underpredicts  the  negative 
peak  pressure.  A  detailed  examination  of  the  Euler  and  Navier-Stokes  solutions  reveals 
strong  directivity.  The  propagated  noise  to  the  sides  of  the  blade  is  relatively  small  and 
non-impulsive.  Upstream  of  the  rotor,  however,  the  propagated  noise  is  very  large  and 
impulsive.  Similar  trends  have  been  observed  in  the  experiments.  Further  examination  of 
the  pressure  contours  in  the  plane  of  the  rotor  shows  the  growth  and  decay  of  the  super¬ 
sonic  pocket  on  the  rotor  blade  surface,  the  initial  separation  of  the  acoustic  wave  from  the 
aerodynamic  field  and  the  propagation  of  the  acoustic  wave  to  the  far-field.  It  is  currently 
impossible  to  obtain  such  detailed  information  from  experiments. 

5.  CONCLUSIONS 

A  free- wake  Euler  and  Navier-Stokes  CFD  method,  called  TURNS,  has  been  developed 
for  helicopter  applications.  This  implicit,  completely  upwind,  finite-difference  numerical 
procedure  of  structured-grids  is  applied  to  calculate  the  viscous  flowfields  of  helicopter  in 
hover,  forward  flight,  and  blade-vortex  interaction  (BVI)  flowfield  of  an  advancing  rotor 
interacting  with  a  concentrated  line  vortex.  The  hovering  flowfields  of  two-bladed  and 
four-bladed  rotors  are  calculated  using  cylindrical  C-H  grid  topology  and  body-fixed  co¬ 
ordinates.  The  vortex  wake  and  its  induced  effects  axe  captured  as  a  part  of  the  overall 
numerical  solution  without  specifying  any  wake  structure  or  position;  that  is,  without  any 
wake  modeling.  The  use  of  periodic  boundary  conditions  with  a  single  blade  enables  the 
construction  of  a  full  helicopter  rotor  flowfield  and  thus  saves  computational  time.  The  cap¬ 
tured  vortical  structure  is  smeared  due  to  grid  coarseness.  Nevertheless,  the  induced  flow 
in  the  plane  of  the  rotor  and  the  near-field  wake  trajectory  axe  computed  fairly  accurately, 
and  the  calculated  surface  pressures  on  the  blade,  thrust,  power,  and  Figure  of  Merit  are 
in  good  agreement  with  experiments.  The  wake  trajectory  is  significantly  improved  when 
new  faxfield  boundary  conditions  axe  introduced. 

Of  the  two  flow  conditions  considered  for  the  two  four-bladed  rotors,  the  BERP  rotor 
produces  more  thrust,  and  the  calculated  Figure  of  Merit  remains  approximately  constant 
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as  Ct  increases.  Results  for  the  UH-60A  and  BERP  blades  at  transonic  high-thrust  con¬ 
ditions  indicate  that  the  BERP  blade  produces  milder  shock-induced  separation  than  the 
UH-60A  blade.  The  BERP  blade  also  shows  a  more  tightly  braided  tip  vortex  structure 
compared  to  the  UH-60A  blade.  Further  improvements  are  needed,  but  the  good  compar¬ 
ison  of  the  UH-60A  rotor  performance  calculations  with  experiments  indicates  that  this 
free-wake  Navier-Stokes  CFD  capability  is  promising  as  a  tool  for  analyzing  the  properties 
of  new  and  exotic  blade  shapes  whose  properties  are  not  known  a  priori.  The  numerical 
method  is  fairly  efficient  and  runs  at  145  MFLOPS  on  the  Cray-YMP  supercomputer.  The 
quasi-steady  Navier-Stokes  calculations  presented  here  for  coarse  and  fine  grids  took  ap¬ 
proximately  1  hour  and  15  hours  of  CPU  time,  respectively,  on  the  Cray-2  supercomputer. 

Overall,  the  study  demonstrates  the  capabilities  of  TURNS  code  in  calculating  the  he¬ 
licopter  rotor  aerodynamic  flowfields,  including  the  acoustics  (high  speed  impulsive  noise), 
in  both  hover  and  forward  flight.  The  aerodynamics  and  acoustics  informations  can  be 
obtained  in  one  single  calculation.  The  agreement  with  experiments  is  very  encouraging, 
demonstrating  the  ability  of  the  solution  scheme  to  capture  the  flowfield  and  acoustic  details 
that  axe  hard  to  obtain  from  experiments. 
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PART  II 


6.  DYNAMIC  STALL  OF  AN  OSCILLATING  WING 
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6.  INTRODUCTION 


The  term  “dynamic  stall”  is  often  referred  to  the  unsteady  separation  and  stall  phe¬ 
nomena  of  aerodynamic  bodies  or  lifting  surfaces  that  are  forced  to  execute  time- dependent 
(unsteady)  motion,  oscillatory  or  otherwise.  It  is  a  complex  fluid  dynamic  phenomenon  of 
practical  importance  and  occurs  on  manuvering  flight  vehicles,  retreating  blades  of  heli¬ 
copter  rotors,  wind  turbine  blades,  and  compressor  cascades.  The  dynamic  stall  phenom¬ 
ena  often  leads  to  the  initiation  of  stall  flutter.  As  summarized  in  extensive  reviews  by 
McCroskey53,54  and  Carr55,  the  majority  of  the  work  on  this  fundamental  fluid  dynamic 
problem  is  devoted  to  the  case  of  airfoils  oscillating  with  moderate  amplitude  in  a  uniform 
freestream.  When  the  airfoil  reaches  fairly  high  angles  of  attack  during  the  oscillatory  cy¬ 
cle,  past  the  static  stall  angle  limit,  the  generated  unsteady  flowfield  is  characterized  by 
massive  unsteady  separation  and  large-scale  vortical  structures.  One  important  difference 
between  this  flowfield  structure  and  that  generated  by  the  static  stall  is  the  large  hystere¬ 
sis  in  the  unsteady  separation  and  reattachment.  The  maximum  values  of  lift,  drag,  and 
pitching-moment  coefficients  can  greatly  exceed  their  static  counterparts,  and  not  even  the 
qualitative  behavior  of  these  can  be  reproduced  by  neglecting  the  unsteady  motion  of  the 
body  surface  (airfoil  or  wing). 

One  of  the  reasons  why  the  flowfield  associated  with  dynamic  stall  is  more  difficult  to 
analyze  than  the  static  stall  is  of  its  dependence  on  a  much  larger  number  of  parameters. 
The  important  ones  are  the  airfoil  shape,  Mach  number,  reduced  frequency  or  pitch  rate, 
amplitude  of  oscillations,  type  of  motion  (ramp  or  oscillatory),  Reynolds  number,  three- 
dimensional  effects,  and  wind  tunnel  effects.  To  date  most  of  the  research  in  this  area  has 
been  performed  for  the  simpler  model  problems  of  two-dimensional  oscillating  airfoils.  Most 
of  what  is  understood  about  the  characteristics  and  various  regimes  of  dynamic  stall  has 
essentially  come  from  experimental  observations,  which  are  mostly  two-dimensional.  At¬ 
tempts  to  calculate  the  quantitative  effects  of  dynamic  stall  abound  in  the  literature  (Refs. 
55-63).  The  purely  laminar  case  is  assumed  to  have  been  solved  (Refs.  57-58),  although 
recent  studies64  show  small-scale  details  of  possible  importance.  However,  the  laminar  cal¬ 
culations  have  not  been  validated  with  experiments  (for  the  lack  of  availability  of  purely 
laminar  data).  The  flows  with  turbulent  boundary  layer  have  not  yet  been  successfully 
solved. 

The  weak  link  in  the  computational/theoretical  methods  for  an  accurate  simulation  of 
these  unsteady  flowfields  is  the  turbulence  modeling.  Of  course,  the  transitional  nature  of 
the  boundary  layer  is  almost  always  neglected;  instead  the  flow  is  approximated  to  be  either 
completely  laminar  or  completely  turbulent  on  the  airfoil  surface.  Such  an  approximation 
may  not  be  correct  if  the  flowfield  is  dominated  by  leading-edge  separation.  In  any  case,  a 
reasonably  good  turbulence  model  must  be  used  to  accurately  calculate  the  nonequilibrium 
nature  of  the  separated  turbulent  boundary  layer  and  the  associated  unsteady  time-lag 
features.  Simple  eddy  viscosity  models,  like  the  Baldwin-Lomax  model26,  have  been  found 
to  be  inadequate.  The  objective  of  this  investigation  is  to  identify  a  reasonably  accurate  and 
robust  turbulence  model.  Severed  turbulence  models  that  are  in  use  in  most  Computational 
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Fluid  Dynamics  (CFD)  codes  are  considered  for  evaluation.  It  should  be  noted  that  Refs.  59- 
62  have  considered  a  similar  exercise. 

The  five  turbulence  models  considered  are  the  Baldwin- Lomax  (B-L)  algebraic  model26, 
the  Renormalization  Group  (RNG)  based  algebraic  model65,  the  half-equation  Johnson- 
King  (J-K)  model66,67,  the  one-equation  Baldwin- Barth  (B-B)  model68,  and  the  one- 
equation  Spalart-Allmaras  (S-A)  model69.  The  performance  of  these  models  is  evaluated 
for  accuracy  and  robustness  by  using  them  to  calculate  the  unsteady,  two-dimensional,  vis¬ 
cous,  flowfields  of  an  oscillating  NACA  0015  wing.  The  accuracy  of  the  calculated  results  is 
determined  by  comparison  with  the  oscillating  wing  experimental  data70  measured  at  the 
U.  S.  Army  Aeroflightdynamics  Directorate  at  NASA  Ames  Research  Center.  The  eventual 
objective  is  to  use  the  turbulence  model  that  calculates  the  unsteady  boundary  layer  and 
flow  physics  accurately  to  simulate  the  three-dimensional  dynamic  stall  of  oscillating  wing 
and  the  retreating  blade  stall  of  a  helicopter  rotor  blade. 

7.  GOVERNING  EQUATIONS  AND  NUMERICAL  SCHEME 

The  governing  equations  considered  are  the  Reynolds-averaged,  two-dimensional, 
Navier-Stokes  equations  in  strong  conservation-law  form.  These  can  be  written  in  a  gener¬ 
alized  body-conforming  curvilinear  coordinate  system  (£,  r?,  r)  as  follows71  : 

drQ  +  d<E  +  a,F  =  ^(d*R  +  d„S)  (15) 


where  £  =  £(x,  y,  t),  rj  =  rj(x,y,t),  and  r  =  t.  Here  (x,  y,  t)  is  the  body-fixed  coordinate 
system;  x  is  along  the  chord  and  y  is  normal  to  the  chord.  Also,  Q  is  the  vector  of  conserved 
flow  variables;  E  and  F  are  the  inviscid  flux  vectors.  These  are  given  by 
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I  (  fmV  +  n,p  \ 

J  I  pvV  +  T)yP  I 
\(e  +  p)V  -  T]tp/ 


The  vectors  R  and  S  axe  the  viscous  stress  vectors  in  the  £  and  tj  directions,  respectively. 
The  viscous  terms  are  retained  in  both  directions  to  resolve  massive  separation  and  these 
are  considered  in  the  thin  layer  approximation.  For  example,  the  vector  in  the  q-direction 
is  written  as 
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where  U  and  V  are  the  contravariant  velocity  components  and  £*,  £v,  tjx,  r)y,  £t,  and  rjt  are 
the  coordinate  transformation  matrics71. 

In  the  above  equations  all  geometrical  dimensions  are  normalized  with  the  airfoil  chord 
length,  c,  the  Cartesian  velocity  components,  u  and  t;  are  scaled  by  the  freestream  sound 
speed  doo,  and  the  time  t  is  normalized  as  tc/a,*,;  P  is  the  static  pressure  normalized  by 
7000;  p  is  the  density  normalized  by  free-stream  density  /><*>;  e  is  the  total  energy  per  unit 
volume  normalized  with  Poo^loi  a  is  speed  of  sound;  Re  is  the  Reynolds  number;  Pr 
is  the  Prandtl  number;  7  is  the  ratio  of  specific  heats;  and  p  is  the  viscosity  coefficient 
normalized  by  its  free-stream  value.  The  pressure  is  related  to  the  density  and  total  energy 
through  the  equation  of  state  for  an  ideal  gas, 

P  =  (7  ~  1)  [e  -  p(u2  +  v2)/2]  (16) 

The  governing  equations  are  written  in  conservation  law  form  for  an  inertial  reference 
frame  ( x0,y0,t ).  Let  r0(<),  v0(t),  and  Q(t)  be  the  position  vector  of  the  origin,  the  velocity, 
and  the  angular  velocity  of  a  non-inertial  frame  (1,  y)  relative  to  the  inertial  frame  (x0,y0). 
The  velocity  v  of  the  non-inertial  frame  relative  to  the  inertial  frame  is: 


v  =  v0(t)  +  Q(t)  x  [r  -  r0(t)j  =  ?„(<) 


When  the  motion  is  reduced  to  pure  rotation  then  v„(t)  =  0.  Positioning  the  origin  of 
the  non-inertial  frame  at  the  origin  of  the  inertial  frame,  for  rotational  type  of  motion, 
makes  r„(t)  =  0.  Therefore,  for  this  restricted  class  of  motion  the  velocity  v  is  expressed 
as,  v  =  Q(t)  x  r(t).  Further  simplification  is  obtained  when  the  body  rotates  around  only 
one  coordinate  axis.  For  example,  when  the  body  rotates  about  the  y-axis  Q(t)  =  (0,u>9) 
then  v  =  u>yex  —  u>yey.  Here  ex  and  ey  are  the  unit  normals  along  x—  and  y-directions, 
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respectively.  The  rotational  speed  is  obtained  from  the  type  of  motion  prescribed  as 
=  da/dt.  The  reduced  frequency  parameter  is  defined  as  either  k  =  de/21/oo,  or 
nfc/Uoo  where  /  is  the  frequency  of  oscillation  and  Uoo  is  the  free-stream  velocity.  Then 
=  ai(2kUoo/c)cos[(2kUoo/c)t]  for  o(<)  =  a0  +  a\sin(Kt)  where  ao  is  the  mean  angle  of 
oscillation,  is  the  amplitude  of  the  pitch,  and  K  =  2kU<x,/c. 

The  numerical  scheme  uses  a  factored,  finite-difference  Beam-Warming  implicit 
algorithm71.  The  viscous  terms  are  retained  in  both  £-  and  ^-directions  and  in  a  thin 
layer  approximation  form.  The  viscous  terms  are  treated  implicitly.  The  approximately 
factored  algorithm  is  given  by 

{/  +  fc[<(Ar„  -  +  (£>,„,),]  }\ 

{<  +  *[<,8?.,  -  Rc-'S, Nfj  +  (!>*.,),] 

(17) 

=  -  Ql,  +  [*<£",  +  «,*£] 

Here  6  is  the  central  difference  operator  and  h  is  the  time-step  that  determines  whether 
the  algorithm  is  first-  or  second-order  time  accurate  .  The  time  index  is  denoted  by  n 
and  AQ'ij  =  (Qfj1  ~  Q?,j)-  The  explicit  inviscid  fluxes  are  given  by  Ejj,  Fij  and  R,j, 

Sij  are  the  explicit  viscous  fluxes.  The  quantities  A ,  B,  M  and  N  in  Eq.  (17)  are  flux 
Jacobian  matrices  obtained  from  the  linearization  of  the  left  hand  side.  Dtmp  and  teD  are 
the  implicit  and  explicit  dissipation  terms.  A  Jameson-type72,  blended  second  and  fourth 
order,  numerical  dissipation  based  on  the  computed  pressure  field  is  used  to  suppress  high 
frequency  numerical  oscillations.  For  subsonic  shock  free  solutions  only  the  fourth-order 
dissipation  is  used,  while  for  transonic  solutions  the  second-order  dissipation  is  activated 
in  the  vicinity  of  shocks  where  the  pressure  jump  has  steep  gradients.  Both  the  implicit 
and  explicit  dissipation  are  scaled  by  the  spectral  radius.  For  the  accuracy  of  calculated 
solutions,  the  added  dissipation  coefficients  are  kept  as  small  as  possible. 

The  errors  introduced  by  the  linearization  and  approximate  factorization  of  the  left 
hand  side  of  the  numerical  algorithm  may  be  minimized  by  performing  Newton  subiterations 
at  each  time-step  during  the  unsteady  calculations.  The  approximation  to  Qn+1  at  each 
subiteration  is  the  quantity  Qp.  When  p  >  2,  during  a  given  level  of  subiteration  to 
convergence,  Qp  =  Qn+1,  but  when  p  =  1  and  no  subiterations  are  performed,  then  Qp  = 
Qn,  and  Qp+1  =  Qn+1 .  In  the  present  study,  the  numerical  experiments  have  demonstrated 
that  because  of  the  small  time-steps  used,  the  Newton  subiterations  are  not  required.  It 
was  also  found  that  the  two  time-level  numerical  scheme  does  not  increase  the  accuracy  of 
the  unsteady  calculations. 

In  the  normal  practice  of  the  thin  layer  approximation  for  viscous  terms,  only  the 
terms  in  the  normal  direction  (5)  sire  retained  because  of  the  large  flow  gradients  in  that 
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direction.  Retaining  the  viscous  terms  for  both  the  directions  (/?,  5),  in  the  present  study, 
was  found  to  be  slightly  beneficial  for  the  deep-stall  cases,  perhaps  because  of  massive 
separation.  For  the  light-stall  case,  however,  the  calculations  performed  retaing  the  viscous 
terms  for  the  77  direction  alone  and  terms  for  both  the  £  and  77  directions  showed  very 
little  difference  between  the  solutions.  Therefore,  the  light-stall  calculations  are  performed 
by  retaining  the  viscous  terms  only  in  the  77  direction.  The  computational  cost  was  not 
significantly  increased  for  doing  this  for  both  the  directions.  Also,  numerical  experiments 
have  demonstrated  that  implicit  treatment  of  the  streamwise  viscous  term,  M  of  Eq.  (17), 
does  not  contribute  to  the  accuracy  of  the  solution  but  results  in  increased  computational 
cost.  Therefore,  implicit  treatment  of  the  streamwise  viscous  term  is  not  used  in  the  present 
study. 

Body-fitted  C-type  computational  grids  are  used  in  all  calculations.  These  are  gen¬ 
erated  using  a  hyperbolic  grid  generator.  The  grids  axe  clustered  at  the  body  surface  in 
the  normal  direction,  leading  edge,  and  trailing  edge  regions.  The  spacing  of  the  first  grid 
point  at  the  surface  in  the  normal  direction  is  0.00002  chord  and  the  grid  boundaries  are 
located  at  15  chords  in  all  directions.  Although  most  of  the  results  presented  are  generated 
using  one  grid  of  size  361x71  with  271  points  on  the  airfoil  in  the  chordwise  direction  and 
the  remaing  in  the  wake  region,  three  other  grids  of  size  181x71,  671x71,  and  361x141  have 
also  been  used  to  study  the  effect  of  grid  size  on  the  solution  accuracy. 

The  orientation  of  the  non-inertial  frame  with  respect  to  the  inertial  frame  is  changing 
at  each  instant  of  time.  Therefore,  after  each  time-step  the  grid  is  moved  to  the  new  location 
and  all  metrics  are  recomputed. 

Boundary  conditions  are  updated  explicitly.  All  flows  are  computed  at  one  subsonic 
free  stream  speed.  For  subsonic  inflow-outflow,  the  flow  variables  at  the  boundaries  are 
evaluated  using  one- dimensional  Riemann  invariant  extrapolation.  At  the  inflow  boundary 
there  is  one  incoming  and  three  outgoing  characteristics.  Therefore,  three  variables,  the 
density,  the  normal  velocity,  and  the  pressure  are  specified  and  the  fourth  variable,  the  axial 
velocity  is  extrapolated  from  the  interior.  At  the  outflow  boundary  there  are  one  incoming 
and  three  outgoing  characteristics  and  only  one  quantity,  the  pressure,  is  specified  while  the 
others  are  extrapolated  from  the  interior.  For  the  density,  a  simple  first-order  extrapolation 
is  used.  On  the  body  surface  a  nonslip  condition  is  applied  for  the  velocities,  viz.,  the 
contravariant  velocities  in  the  body-fixed  coordinates  are  set  to  zero.  It  should  be  noted 
that  the  surface  velocity  is  non-zero  because  of  the  body  motion  through  the  unsteady 
metrics.  The  boundary  layer  approximation  is  used  to  obtain  the  surface  pressures  as 
Pi,  1  =  Pi, 2-  For  C-type  grids  used  in  this  study  averaging  of  the  flow  variables  at  the  wake 
cut  is  used. 
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8.  TURBULENCE  MODELS 


All  flows  computed  in  this  study  are  assumed  to  have  fully  turbulent  boundary  layer 
on  both  the  upper  and  lower  surface  of  the  airfoil  by  neglecting  the  laminar  and  trasitional 
boundary  layer.  In  the  experiments  of  Piziali70,  with  which  all  present  calculations  are 
compared,  the  boundary  layer  is  tripped  in  the  leading-edge  region.  For  turbulent  viscous 
flows,  the  non-dimensional  viscosity  fi  in  the  viscous  flux  vectors  is  calculated  as  the  sum 
total  of  the  laminar  and  turbulent  viscosity.  It  is  the  determination  of  this  turbulent  vis¬ 
cosity  that  is  of  special  significance  and  the  focal  point  of  this  present  study.  As  mentioned 
before,  five  different  turbulence  models  are  used  for  calculating  turbulent  eddy  viscosity 
and  the  unsteady  flowfield  of  an  NACA  0015  oscillating  wing.  The  results  are  used  to 
evaluate  their  performance.  The  details  of  how  these  models  sure  developed  are  described 
elsewhere26’65-69.  The  following  paragraphs  describe  briefly  the  salient  features  of  these 
models  and  the  specific  versions  used  in  the  present  investigation. 


8.1  Baldwin- Lomax  (B-L)  model 


This  is  a  two-layer,  zero-equation  model.  It  is  patterned  after  Cebeci-Smith  model73 
and  introduces  a  modification  that  eliminates  the  need  to  search  for  the  edge  of  the  bound¬ 
ary  layer  to  determine  length  scale.  It  is  the  most  commonly  used  turbulence  model  avail¬ 
able  in  most  of  the  CFD  codes.  Its  strength  and  weakness  are  well  known  in  CFD  commu¬ 
nity;  it  predicts  accurately  the  steady  flows  with  little  or  no  separation  and  fares  poorly  if 
there  is  large  separation,  either  shock-induced  or  otherwise.  It  uses  an  inner  and  outer  layer 
formulation  for  determining  the  turbulent  viscosity  with  a  smooth  transition  that  spreads 
over  the  two  regions.  It  uses  a  classical  mixing- length  hypothesis  for  the  inner  layer  with  a 
van  Driest  damping  function  to  force  the  eddy  viscosity  at  the  wall  to  zero.  In  the  outer 
layer,  the  length  scale  is  fixed  by  the  location  where  the  product  of  distance  from  the  solid 
wall  and  vorticity  reaches  a  maximum  in  the  boundary  layer.  The  Klebanoff’s  intermit- 
tency  factor  is  used  to  drive  the  eddy  viscosity  to  zero  in  the  outer  flow  away  from  the  wall. 
Some  of  the  constants  of  the  theory  Eire  determined  by  correlating  with  experimental  data. 
The  details  of  the  theory  are  described  in  Ref.  26. 


8.2  RNG  model 


Another  algebraic  eddy  viscosity  model  was  proposed  recently,  for  the  closure  of  the 
Reynolds  averaged  Navier-Stokes  equations,  based  on  the  Renormalization  Group  (RNG) 
theory  of  turbulence65.  The  algebraic  model,  although  free  from  uncertainties  related  to 
the  experimental  determination  of  empirical  modeling  constants,  still  requires  specification 
of  an  integral  length-scale  of  turbulence,  similar  to  the  B-L  model,  which  reduces  the 
generality  of  the  model.  In  this  model  the  integral  scale  is  assumed  to  be  proportional  to 
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the  boundary  layer  thickness  6,  and  the  eddy  viscosity  is  obtained  as  in  Ref.  65  from 


v  =  I/| 


1  + 


(18) 


where  v  =  vt  +  i/j,  the  subscripts  t  and  l  refer  to  the  turbulent  and  laminar  components, 
respectively  and  6  is  the  boundary  layer  thickness.  It  is  determined  from  8  =  1.2  yj/2  where 
Di/2  is  the  normal  distance  from  the  wall  at  which  the  vorticity  function  F(y)  (see  Ref.  26) 
attains  its  half  amplitude67.  H  is  the  Heaviside  step  function  and  4>  is  the  dissipation 
rate,  it  is  determined  by  assuming  production  equals  dissipation  for  equilibrium  flows. 
The  parameter  a  =  0.0192  corresponding  to  the  von  Karman  constant  k  =  0.372  and  the 
parameter  Cc  =  75.  The  turbulent  eddy  viscosity  is  then  obtained  by  solving  Eq.  (18)  at 
every  point  in  the  computational  domain.  In  estimating  the  eddy  viscosity  with  this  model 
in  this  study,  the  model  is  applied  only  to  the  suction  side  of  the  airfoil  (upper  surface) 
while  the  pressure  side  (having  attached  flow)  and  wake  regions  are  computed  with  the  B-L 
model.  Application  of  the  model  to  both  the  upper  and  lower  surfaces  essentially  gave  the 
same  results  as  the  one  obtained  by  applying  for  only  uppper  surface;  so  the  latter  was 


used. 


8.3  Johnson-King  (J-K)  model 


The  above  two  models,  viz.,  the  B-L  and  RNG  models,  are  termed  equilibrium  models 
meaning  that  the  eddy  viscosity  instantaneously  adjusts  to  the  local  flow  without  any 
history  effects.  The  next  three  models  presented  are  called  non-equilibrium  models  in 
which  the  calculated  eddy  visocity  accounts  for  the  upstream  history  of  the  flow. 


From  the  time  Johnson  and  King  first  introduced  their  half-equation  turbulence 
model66,  there  have  been  several  modifications  and/or  enhancements  to  improve  their  orig¬ 
inal  model  for  separated  flows61.  In  the  present  application  to  unsteady  flows,  the  original 
version  of  this  model  is  used66  and  is  briefly  described  in  the  following  paragraphs. 

The  Johnson-King  model66  takes  into  account  the  convection  and  diffusion  effects  on 
the  Reynolds  shear  stress  —u'w'  in  the  streamwise  direction.  The  eddy  viscosity  is  given 
by 

Vt  =  vu  [l  -  exp{-^-) j  (19) 

where  vti ,  uto  describe  the  eddy  viscosity  variation  in  the  inner  and  outer  part  of  the  bound¬ 
ary  layer.  The  inner  eddy  viscosity  is  computed  as 

vti  =  D2KyyJ(PJtf)max 
D  =  1  -  e~{*/A+) 

where  the  constant  A+  =  15.  The  outer  eddy  viscosity  is  given  by 
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ut,  =  *(*)  ■  (0.0168f7e6*7) 


(21) 


where  6*  is  the  boundary  layer  displacement  thickness,  7  is  the  KlebanofTs  intermittency 
function  given  by  7  =  [l  +5.5(  J)  ]  1 ,  and  a(x)  is  obtained  from  the  solution  of  an  ordinary 
differential  equation  which  describes  the  development  of  —  u#u>'|maz  along  the  path  of  the 
maximum  shear  stress.  The  effects  of  convection  and  diffusion  on  the  Reynolds  shear  stress 
development  are  accounted  from  the  solution  of  the  following  ordinary  differential  equation 
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dx 
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a, «[0.7 -({)„] 

Here  Cdif  and  aj  are  modeling  constants  um  is  the  maximum  average  mean  velocity  and 


g  =  [— u'tw'n,]  ,  and 

Q* V  =  l— u<u,/)m,e4f] 

where  Lm  is  the  dissipation  length  evaluated  as 


Lm  =0.40y  for  ym/6  <  0.225  . 

Lm  =0.096  for  ym/6  >  0.225  ( 

The  boundary  layer  thickness,  6,  is  determined  in  the  same  way67  as  explained  in  the 
discussion  of  RNG  model.  The  equilibrium  shear  stress  gtq  in  Eq.  (22)  is  determined  from 
the  following  equilibrium  eddy  viscosity  distribution 


^.e,  =*'«..€»[  1  *■)  I 

=D2*yyJ(Lu'w,)m,tq  ^ 

vt.,cq  =0.0168^6*7 

where  Ut  is  the  velocity  at  the  edge  of  the  boundary  layer. 

An  implicit  Euler  method  is  used  for  the  numerical  solution  of  Eq.  (22),  and  the 
maximum  shear  stress  at  each  iteration  level  is  updated  as  follows 


It  should  be  noted  that  the  unsteady  term  is  neglected  in  the  above  formulation.  Solutions 
with  the  Johnson-King  turbulence  model  are  obtained  as  follows.  First  a  convergent  solu¬ 
tion  using  the  Baldwin-Lomax  turbulence  model  for  the  entire  flowfield  is  obtained.  Then 
the  Johnson-King  model  is  applied  only  to  the  upper  surface  of  the  airfoil  as  using  it  for 
both  the  surfaces  did  not  change  the  results.  To  initiate  the  solution  <j(x)  in  Eq.  (21)  is 
set  unity  and  it  is  allowed  to  change  according  to  Eq.  (25).  It  should  be  noted  that  the 
Johnson-King  model  reduces  to  the  Cebeci-Smith  model73  when  <7(1)  is  identically  equal 
to  one. 

8.4  One— Equation  Models 

The  B-L  and  the  RNG  models  are  equilibrium  models,  in  which  the  production  is 
identically  equal  to  the  dissipation.  The  J-K  model  is  an  improvement  over  the  equilib¬ 
rium  turbulence  models  because  it  accounts  for  the  evolution  of  the  maximum  shear  stress 
through  the  solution  of  an  ordinary  differential  equation  (ODE).  It,  therefore,  attempts  to 
calculate  the  non-equilibrium  turbulent  boundary  layer.  The  validity  of  these  models  is  lim¬ 
ited  and  questionable  when  applied  to  a  flow  environment  consisting  of  unsteady  separated 
flow  with  multiple  shear  layers. 

Recently,  several  one-equation  models  have  been  developed  for  use  in  place  of  these 
lower-order  turbulence  models.  In  the  present  investigation  two  such  models  are  consid¬ 
ered  for  investigation.  These  are  the  Baldwin- Barth68  and  Spalart-Allmaras69  models. 
The  primary  advantage  of  these  methods  is  that  they  do  not  require  the  evaluation  of 
flow-dependent  length  scales,  such  as  the  boundary  layer  thickness.  The  validity  of  these 
models  for  steady  flows  has  been  demonstrated,  but  only  in  a  limited  sense.  In  the  present 
investigation  these  models  are  tested  for  several  unsteady  attached  and  separated  flows  over 
oscillating  airfoils. 

8.4.1  Baldwin-Barth  (B-B)  model 

This  one-equation  model68  is  derived  from  the  simplified  form  of  the  k  —  e  model 
equations.  It  solves  for  the  modified  turbulent  Reynolds  number  uRt  from 

={c.,h  -  c„)\] vRtP 

+  (p+  —)V2(„Rt)  (26) 

-  —  (Vi/t)  •  V(uRt) 

This  is  a  partial  differential  equation  for  the  field  quantity  Rt  =  k2  jvt  =  RrfziRr )> 
named  turbulent  Reynolds  number.  The  turbulent  Reynolds  number  is  related  to  the  eddy 
viscosity  as 
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Vt  =  vcpf^Rr  =ucHftlf3RT 


(27) 


where 
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Here  y+  =  ury/t/  and  ur  is  the  skin  friction  velocity.  The  constants  used  for  the  B-B  model 
are  the  same  as  in  their  original  paper68  and  are  given  by: 

k  =0.41,  ctl  =  1.2,  c«,  =  2.0 

Cp  =0.09,  A+  =  26.,  A2+  =  10. 

This  model  is  applied  to  the  entire  flowfield  to  compute  the  eddy  viscosity. 


8.4.2  Spalart-Allmaras  (S-A)  model 

The  second  one-equation  model  used  in  the  present  investigation  is  the  Spalart- 
Allmaras  (S-A)  model69.  This  model  requires  the  solution  of  a  transport  (partial  dif¬ 
ferential)  equation  for  the  turbulent  eddy  viscosity.  This  equation  was  constructed  using 
empirical  criteria  and  arguments  from  dimensional  analysis.  It  has  many  similarities  with 
the  B-B  model  but  it  is  relatively  simpler  for  its  numerical  implementation.  The  S-A  model 
also  incorporates  transition  location  specification,  although  the  present  investigation  treats 
the  boundary  layer  as  turbulent  on  the  entire  surface.  The  eddy  viscosity  is  obtained  from 
the  solution  of  the  following  partial  differential  equation. 

=c» i(l  -  ft2)Su  +  i  [V  •  ((*/  +  i/)V«/)  +  cM(V«/)2 
—  jcwi/u,  —  [~j]  +  fti&u1 


34 


where 


fti  =ctigtexp(  -  +  92M) 

ft2  =ct3exp(-ct4X2)  ^ 

gt  =min(0.1,  AU/u>trAx) 

where  x  =  vfv  and  uiT  is  used  here  to  denote  the  vorticity  at  the  wall  at  the  boundary 
layer  trip  point.  The  constants  of  this  model  have  been  chosen  the  same  as  in  the  original 
reference69,  and  the  transition  location  was  set  at  the  airfoil  leading  edge.  The  model 
constants  are: 


Cil 

=0.1355 

c62 

=0.622 

a 

=2/3 

CU)1 

=cj, i/k2  H 

Cu,2 

=2.0 

Cw3 

=0.3 

K 

=0.41 

—  /7  | 

(  1  +  cw2 

—g\ 

V  +  <i, 

) 
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g  =r  +  cw2(r*  -  r),r  = 


vt 


Sk2cP 


cti  =1.0 

c«2  =2.0 
Ct3  =1-1 

Ct4  =2.0 


The  relative  computational  costs  for  each  of  the  five  turbulence  models  is  discussed 
towards  the  end  of  Results  Section. 


9.  RESULTS  AND  DISCUSSION 


This  paper  describes  an  attempt  to  evaluate  five  different,  but  widely  used,  turbulence 
models  in  calculating  the  unsteady,  two-dimensional  flowfields  of  an  oscillating  NACA  0015 
wing  spanning  a  wind  tunnel  wall  and  an  end  plate.  The  test  cases  considered  for  the 
calculation  correspond  to  the  wind  tunnel  conditions  of  an  experiment  of  a  NACA  0015 
oscillating  wing  conducted  in  the  7x10  foot  wind  tunnel  of  the  U.  S.  Army  Aeroflightdy- 
namics  Directorate,  located  at  NASA  Ames  Research  Center19.  The  flow  conditions  of  the 
experiment  are  as  follows:  free-stream  Mach  Number,  Moo  =  0.29  and  Reynolds  number, 
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Re  =  1.95xl06  based  on  the  chord  of  the  wing  and  free-stream  speed.  Four  mean  angles  of 
<*o  =  4°,  11°,  15°,  and  17°  are  considered  with  an  amplitude  of  pitch  of  aj  =  ±4.2°  around 
the  mean  angle  and  a  reduced  frequency,  k  =  0.1.  The  instantaneous  angle  of  attack  of  the 
wing  is  given  by  a(t)  =  <*o  +  sin(Kt). 


The  present  calculations  are  done  for  free-air  conditions  and,  therefore,  the  wind  tunnel 
walls  are  not  included  in  the  calculations.  The  unsteady  calculations  for  oscillating  airfoils 
are  usually  started  from  the  steady  state  solution  at  the  lowest  angle  of  attack  of  the  pitch 
cycle.  However,  to  check  the  accuracy  of  the  solution  method,  steady-state  solutions  were 
calculated  for  several  angles  of  attack  in  a  time-accurate  manner.  Figure  37  shows  sample 
results  of  steady  surface  pressures  for  a  =  13°  (Fig.  37a)  and  17°  (Fig.  37b)  compared 
with  experiments.  The  B-L  turbulence  model  yielded  satisfactory  surface  pressures  for  the 
case  of  a  =  11°  where  the  flow  is  essentialy  attached  (not  shown  here).  For  the  mildely 
separated  case  of  a  =  13°  and  the  massively  separated  cases  of  a  =  15°  and  17°,  the  B-L 
model  predicted  surface  pressures  that  had  higher  leading-edge  suction  peaks  and  milder 
separation.  It  gave  higher  lift  and  lower  drag  values.  Therefore  more  accurate  turbulence 
models  were  needed  to  predict  satisfactory  airloads.  As  seen  in  Fig.  37,  the  J-K  model 
is  able  to  predict  the  surface  pressures  more  accurately  than  the  B-L  algebraic  model  at 
these  conditions.  Although  not  shown  here,  the  B-B  one-equation  model  fared  as  good  at 
13°  and  slightly  better  at  17°  than  the  J-K  model  in  predicting  the  overall  quasi-steady 
airloads.  The  airfoil  is  completely  stalled  at  the  a  =  17°  condition. 

Table  1  lists  the  calculated  force  coefficients  for  these  two  quasi-steady  flow  conditions 
and  the  experimental  measurements.  The  airloads  are  calculated  in  this  study  by  integrat¬ 
ing  the  surface  pressures.  Therefore,  the  drag  coefficient  shown  is  that  due  to  pressure  drag 
only.  There  are  large  fluctuations  in  airloads  in  the  experimental  data  not  only  from  cycle 
to  cycle  at  a  single  station,  but  also  from  station  to  station  for  a  two-dimensional  wing. 
The  results  listed  for  the  coefficients  of  lift  (Ci),  drag  (Cd),  and  pitching-moment  about 
quarted-chord  (Cm)  in  this  table  are  time- averaged  values  over  a  large  period  of  time  for 
the  calculations.  The  experimental  value  was  read-off  from  Ref.  70.  At  the  lower  angle  of 
a  =  13°,  the  J-K  and  B-B  models  did  better  in  predicting  the  overall  airloads  than  the 
B-L  model,  although  all  of  them  fail  to  predict  accurately  the  drag  and  pitching-moment. 
At  the  stalled  condition  of  17°  the  flow  was  unsteady.  The  results  are  in  poor  agreement 
with  each  other  and  with  experiment,  although  the  agreement  with  experiments  improved 
as  the  models  got  more  sophisticated  from  B-L  to  B-B  models. 


The  discussion  of  unsteady  results  is  divided  into  three  flow  regimes,  viz., 
a)  attached  flow  corresponding  to  <*o  =  4°;  b)  light-stall  case  corresponding  to 
oo  =  11°;  toid  c)  deep-stall  cases  corresponding  ao  =  15°  and  17°.  After  calculating 
the  quasi-steady  solution  at  the  lowest  angle  of  the  pitch  cycle  for  each  of  the  condition, 
the  airfoil  is  made  to  execute  pitching  oscillations  rotating  about  its  quarter-chord  point. 
The  unsteady  flow  evolution  is  monitored.  Most  of  the  results  presented  in  the  following 
paragraphs  were  calculated  using  361x71  grid  with  10,000  constant  time-steps  per  osciallat- 
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ing  cycle.  This  number  of  time-steps  per  cycle  corresponds  to  a  nondimensional  time-step 
of  At  =  0.0108,  based  on  c  and  Uoo-  An  explicit  dissipation  coefficient  of  ee  =  0.05  was  used 
as  the  standard  value.  A  parametric  study  of  the  grid  size,  the  value  of  ce,  and  number  of 
time-steps  per  cycle  has  also  been  done  for  one  turbulence  model,  for  the  light-stall  case, 
to  identify  optimum  values  to  be  used  for  all  calculations  and  this  discussions  is  presented 
in  the  section  on  the  Light-stall. 

9.1  Unsteady  attached  flow  case:  a(t)  =  4°  +4.2 °sin(Kt) 

This  flow  condition  essentially  serves  to  validate  the  accuracy  of  the  flow  solver  in 
calculating  unsteady  attached  flow.  Figure  38  presents  the  calculated  unsteady  airloads 
with  different  turbulence  models  compared  to  the  experimental  data70.  The  experimen¬ 
tal  data  is  averaged  over  20  cycles  period  and  the  calculations  are  not  averaged  but  were 
found  to  be  repeatable  beyond  the  second  cycle.  The  results  of  calculations  with  different 
turbulence  models  compare  favorably  with  each  other  for  lift  and  pressure  drag,  but  the 
pitching-moment  appears  to  be  more  sensitive.  The  RNG,  J-K,  and  S-A  models  give  rel¬ 
atively  good  comparison  among  themselves  and  with  experimental  data  for  lift,  drag,  and 
pitching-moment.  It  is  surprising  to  see  that  the  B-B  model  does  relatively  poorly  com¬ 
pared  to  the  above  three  models  for  lift  and  pitching-moment.  Similarly,  the  B-L  model 
does  poorly  in  predicting  the  lift  and  pithing-moment.  It  consistently  predicts  higher  lift 
and  lower  pitching-moments  compared  to  the  rest  of  the  models.  It  should  be  noted  that 
the  scales  used  in  presenting  the  airloads  in  Fig.  38  are  expanded  to  bring  out  the  differ¬ 
ences  clearly  for  various  turbulence  models,  but  in  the  scales  used  in  rest  of  the  study  for 
different  airloads,  the  results  are  well  within  the  range  of  experimental  scatter.  The  trends 
of  the  calculated  results  show  that  the  if  the  results  are  tilted  slightly,  they  perhaps  agree 
better  with  experiments.  As  observed  by  McCroskey  and  Pucci74  in  their  experiments  with 
NACA  0012  oscillating  wing,  such  a  trend  can  be  attributed  to  wind  tunnel  wall  interference 
effects. 

Figure  39  shows  the  details  of  unsteady  pressure  distributions  by  harmonic  com¬ 
ponents,  where  the  decomposition  of  Cp  is  according  to  Cp  =  CPo  +  CPllsin(Kt)  + 
ccos(Kt)  +  CP7sin(2Kt  +  <fo)  with  4>i  as  the  phase-shift.  In  order  to  stretch  out  the 
leading-edge  region  of  the  airfoil,  the  various  harmonics  axe  plotted  against  y/x.  Such  a 
representation  brings  out  the  variations  and  discrepancies  better  with  linear  theory.  Note 
that  the  differences  seen  in  the  lift  and  pitching-moment  hysteresis  curves  for  B-L  and  B-B 
models  in  Fig.  38  relatively  translate  into  small  differences  in  the  mean  and  quadrature 
components  in  Fig.  39.  No  experimental  pressure  data  was  available  for  comparison  with 
these  results,  but  comparison  of  the  four  components  with  the  measurements  of  Ref.  74  for 
NACA  0012  oscillating  airfoil  for  attached  unsteady  flow  case  show  similar  behavior  for  all 
the  components. 


37 


9.2  Light-stall  case:  a(t)  =  11°  +  4.2°sm(A't) 


Figures  40-45  show  the  results  for  the  light-stall  case  for  ao  =  11°-  In  order  to  deter¬ 
mine  the  optimum  values  of  time-step,  At,  the  explicit  numerical  dissipation  coefficient,  ce, 
and  size  of  the  grid  for  a  reasonably  accurate  solution,  parametric  calculations  have  been 
done  using  one  turbulence  model  (J-K)  and  varying  these  variables  one  at  a  time  while 
keeping  the  others  constant.  For  example,  Fig.  40  presents  results  of  unsteady  airloads  hys¬ 
teresis  for  three  values  of  At  for  a  grid  size  of  361x71  and  ee  =  0.05  along  with  experimental 
data.  In  Fig.  40,  At  =  0.0108  corresponds  to  using  10,000  time-steps  per  cycle.  Similarly, 
At  =  0.0216  corresponds  to  5,000  time-steps  per  cycle.  The  maximum  number  used  in 
the  present  calculations  is  20,000  time-steps  per  cycle  for  deep-stall  cases.  This  contrasts 
with  50,000  to  100,000  time-steps  per  cycle  used  by  other  investigators62,63  in  similar  situa¬ 
tions.  The  large  number  of  time-steps  in  these  investigations  was  dictated  by  the  numerical 
stability  of  their  solution  method.  Obviously  these  become  prohibitively  expensive  when 
extended  to  three-dimensions.  As  seen  in  Fig.  40,  a  A t  =  0.0108  corresponding  to  10,000 
time-steps  per  cycle  is  a  good  compromise,  and  this  is  the  number  used  for  calculating  all 
the  results  presented  in  this  study. 

The  results  of  unsteady  airloads  presented  in  Fig.  41  shows  that  all  three  forces  are 
very  sensitive  to  the  dissipation  coefficeint  ce.  These  airloads  are  calculated  using  a  A t  = 
0.0108  for  a  361x71  grid  with  J-K  turbulence  model.  For  the  four  values  of  ee  used,  the 
upstroke  results  of  lift,  drag,  and  pitching- moment  are  nearly  the  same.  Differences  are  seen 
for  the  downstroke  results  with  higher  values  of  ce  giving  better  drag  and  pitching-moment 
predictions.  From  this,  a  value  of  ce  =  0.05  is  chosen  as  the  reference  value  for  the  rest 
of  the  calculations.  For  a  value  less  than  this,  the  solutions  at  other  higher  mean  angles 
showed  oscillatory  behavior  for  the  unsteady  airloads. 

Figure  42  presents  the  unsteady  airloads  results  calculated  on  different  grid  sizes. 
Again,  J-K  model  is  used  along  with  ce  =  0.05  and  A t  =  0.0108  for  these  calculations.  As 
seen  the  solutions  are  quite  sensitive  to  the  grids  used.  It  should  be  noted  that  all  grids 
have  the  same  wall  spacing  of  the  first  grid  point  in  the  normal  direction  and  the  boundaries 
are  located  at  the  same  distance  from  the  airfoil.  It  is  surprising  to  see  that  even  fine  grids 
have  produced  poor  results  compared  to  the  361x71  grid.  This  perhaps  is  due  to  using  the 
same  €e  for  all  cases.  It  is  evident  from  the  results  presented  in  Fig.  42  that  for  the  ee  and 
A t  used,  the  361x71  grid  appears  to  be  the  best  choice  and  this  is  the  grid  used  for  most 
results  presented  in  this  study. 

A  similar  result  of  grid-sensitive  study  using  B-B  model  is  presented  in  Fig.  43.  For 
this  turbulence  model,  the  unsteady  airloads  show  less  sensitivity  to  different  grids.  These 
results  are  also  calculated  using  the  same  values  of  At  and  ce  as  those  of  Fig.  42  for  J-K 
model.  As  seen,  the  drag  and  pitching-moment  have  better  agreement  with  experiment  for 
both  upstroke  and  downstroke,  but  the  lift  on  the  downstroke  has  very  poor  agreement  with 
experiment  indicating  that  the  flow  reattachment  is  not  complete  until  after  the  upstroke 
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begins.  Again  the  361x71  grid  appears  to  be  the  optimum  grid  for  the  given  At  and 
dissipation. 


Using  the  above  suggested  optimum  values  for  At ,  ee,  and  316x71  grid,  unsteady 
flowfield  solutions  are  calculated  using  all  the  turbulence  models.  Figure  44  presents  the 
unsteady  airloads  results  from  these  solutions.  Comparison  of  results  for  different  turbu¬ 
lence  models  reveal  that  every  model  behaves  differently.  The  chief  characteristic  of  all 
these  solutions  is  that  all  models  produce  trailing-edge  separation.  The  B-L  model  (not 
shown  in  this  figure)  produces  the  least  separation.  The  lift,  drag,  and  pitching-moment 
hysteresis  curves  for  this  model  are  distinctly  different  from  the  rest  of  the  solutions  and, 
in  particular,  the  lift  and  pitching-moment  are  in  poor  agreement  with  experiment.  The 
lift  hysteresis  curves  for  the  J-K,  RNG,  and  S-A  models  are  in  good  agreement  with  exper¬ 
iments,  although  the  RNG  model  has  slightly  higher  lift  during  the  upstroke.  In  general, 
all  models  show  good  agreement  with  each  other  and  with  experiment  for  the  unsteady  air¬ 
loads  on  the  upstroke  part  of  the  cycle.  As  seen  here,  no  one  model  has  perfect  agreement 
with  experiment  for  all  airloads.  Once  the  flow  separates  on  the  upper  surface  of  the  airfoil, 
before  the  downstroke  begins,  the  B-B  model  appears  to  recover  very  slowly  to  the  attached 
flow  condition  compared  to  the  other  models.  In  fact  the  boundary  layer  reattachment  is 
not  complete  for  the  B-B  model  until  after  the  upstroke  begins.  Therefore,  the  lift  stays 
very  low  until  the  upstroke  begins.  Although  not  apparent  in  the  drag  curve,  this  is  very 
well  seen  in  the  pitching-moment  curve. 

The  RNG  and  J-K  models  produce  nearly  the  same  extent  of  separation  that  is  much 
larger  than  what  the  B-B  model  produces.  They  also  appear  to  have  similar  recovery  in 
the  downstroke  much  better  than  the  B-B  model.  Both  models  predict  very  similar  drag 
and  pitching-moment  that  are  in  poor  agreement  with  experiment  for  the  downstroke.  The 
S-A  model  produces  separation  similar  in  extent  to  the  B-B  model,  but  has  good  recovery 
just  like  the  J-K  model.  The  lift  hysteresis  for  this  model  is  in  good  agreement  with 
experiment,  like  the  J-K  model,  and  the  drag  and  pitching-moment  are  in  better  agreement 
with  experiment  than  the  J-K  model  but  not  as  good  as  the  B-B  model.  All  the  models 
except  the  B-B  model  predict  poor  pitching-moments.  The  B-B  model  produces  nearly  the 
right  amount  of  separation  and,  that  is  the  reason,  it  predicts  pitching-moment  correctly 
including  its  cusp-like  behavior  at  the  end  of  upstroke.  The  hump-like  behavior  at  the  end 
of  downstroke  is  the  result  of  poor  recovery  as  discussed  before. 


In  general,  all  models  predict  the  unsteady  airloads  reasonably  well  for  the  upstroke 
and  they  all  behave  differently  during  the  downstroke.  The  B-B  model  is  the  only  model 
that  predicts  the  pitching-moment  correctly  and  in  agreement  with  experimental  data  ex¬ 
cept  for  the  part  at  the  end  of  downstroke.  But  it  nicely  recovers  as  the  upstroke  begins. 
No  leading-edge  separation  is  seen  for  any  of  the  models,  perhaps  due  to  negleting  laminar 
and  transition  regions  and  assuming  turbulent  boundary  layer  for  the  entire  airfoil. 

Figure  45  shows  the  harmonic  components  of  unsteady  pressures  calculated  for  this 
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light-stall  case  using  the  J-K  and  B-B  models.  All  four  parts  of  this  figure  axe  different 
from  those  of  the  preceding  attached  case,  q0  =  4°,  shown  in  Fig.  39.  The  large  changes 
seen  in  these  curves  can  be  attributed  to  the  nonlinear  behavior  of  the  flow  at  this  mean 
angle  of  ao  =  11°.  The  two  models  predict  very  similar  mean  and  in-phase  components. 
But  the  very  different  pitching-moments  produced  by  the  two  models  is  apparent  in  the 
out-of-phase  component. 

9.3  Deep-stall  case:  a(t)  —  ao  +  4.2°  sin(Kt) 

In  contrast  to  the  light-stall  with  mild  trailing-edge  separation  seen  for  ao  =  11°  case, 
the  deep-stall  cases  for  ao  =  15°  and  17°  are  dominated  by  massive  flow  separation  and 
highly  nonlinear  flow  behavior.  The  separation  that  originates  in  the  trailing-edge  region, 
during  the  upstroke  of  the  oscillatory  cycle,  continues  to  spread  upstream  as  the  airfoil 
motion  changes  to  downstroke.  The  unsteady  stall  behavior  in  this  regime  is  characterized 
by  the  shedding  of  a  large  vortex-like  structure  during  the  downstroke  of  the  cycle.  This 
structure  convects  over  the  upper  surface  of  the  airfoil  and  leaves  the  trailing  edge  before 
the  downstroke  part  of  the  oscillatory  cycle  is  completed.  As  a  result,  unsteady  airloads 
far  in  excess  of  the  static  counterpart  are  generated  during  the  upstroke  and  large  amounts 
of  hysteresis  occur  during  the  remainder  of  the  cycle.  The  scale  of  the  viscous-inviscid 
interaction  zone  is  also  large,  producing  a  viscous  layer  thickness  of  the  order  of  the  airfoil 
chord,  particularly  during  the  vortex  shedding  process. 

Figure  46  shows  the  calculated  lift,  drag,  and  pitching- moment  hysteresis  loops  for 
the  ao  =  15°  case  using  different  turbulence  models.  The  results  are  also  compared  to  ex¬ 
perimental  data.  For  reference,  the  B-L  results  are  also  shown  for  this  case.  Examination 
of  individual  curves  reveals  the  following  behavior.  The  lift  hysteresis  is  in  general  agree¬ 
ment  with  calculations  for  all  turbulence  models  except  for  the  B-L  model  which  shows 
an  oscillatory- type  behavior  for  *he  loads  during  downstroke.  All  models  show  very  good 
agreement  with  each  other  and  w!„h  the  experiment  for  the  upstroke  lift  curve.  There  are 
differences  in  the  size  of  the  dynamic  stall  vortex  produced  by  each  model  and  its  convection 
downstream.  Nevertheless,  there  are  only  minor  differences  in  the  downstroke  lift  curves 
for  these  models.  As  before,  the  B-B  model  shows  clearly  that  it  is  slow  in  the  recovery  and 
reattachment  process.  As  a  result,  the  lift  stays  lower  compared  to  the  experimental  value 
and  the  recovery  is  not  complete  until  the  airfoil  changes  its  attitude  from  downstroke  to 
upstroke  motion. 


Reasonbly  good  agreement  of  lift  curves  with  experiment  is  not  an  indication  to  how  the 
drag  and  pitching-moment  are  predicted.  Unlike  the  light-stall  results  for  drag  and  pitching- 
moment  of  Fig.  44,  the  deep-stall  case  shows  a  steep  rise  in  drag  curve  and  negatively  steep 
rise  in  the  pitching- moment  towards  the  end  of  the  upstroke  of  the  cycle.  The  J-K  model 
which  predicted  good  lift  and  drag  curves  for  mean  angles  ao  =  4°  and  11°,  also  shows 
good  agreement  for  lift  hysteresis  with  experiment,  but  predicts  poorly  for  the  drag  and 
pitching-moment.  Examination  of  instantaneous  particle  flow  pictures  of  flow  at  different 


40 


airfoil  positions  during  the  downstroke  reveals  that  the  dynamic  stall  vortex  leaves  the 
surface  much  sooner  than  what  other  models  predict.  Except  for  the  B-L  and  J-K  models, 
all  other  models  have  good  qualitative  agreement  with  experiment  for  drag  and  pitching- 
moment. 

The  behavior  of  the  intantaneous  surface  pressures  during  the  oscillatory  cycle  is  shown 
in  Fig.  47  for  the  upper  surface  of  the  airfoil  for  the  B-B  model.  As  seen  in  this  figure,  the 
leading-edge  suction  peak  continues  to  rise  through  the  upstroke  without  stall.  The  peak 
angle  of  attack  at  the  end  of  upstroke  is  19.2°.  This  angle  is  about  6°  above  the  static 
stall  angle  for  this  airfoil,  as  unsteady  effects  extend  the  dynamic  lift  beyond  the  static 
stall  angle.  The  leading-edge  suction  peak  suddenly  collapses  immediately  after  the  airfoil 
starts  the  downstroke,  as  revealed  by  the  surface  pressure  distributions.  Another  revealing 
feature  of  this  plot  is  that  the  vortex  shedding  phenomenon  manifests  itself  in  the  pressure 
distributions  on  the  downstroke. 

The  harmonic  components  of  the  unsteady  pressures  for  the  B-B  model  are  presented 
in  Fig.  48.  Three  of  the  four  components  are  very  different  from  those  of  the  light-stall 
case  presented  in  Fig.  45  and  also  of  the  unseparated  flow  presented  in  Fig.  39.  Only  the 
mean- component  is  qualitatively  similar  in  shape  to  Fig.  45a.  The  in-phase,  out-of-phase, 
and  the  second  harmonic  are  all  changed  by  the  massive  separation  and  the  presence  of  a 
large-scale  dynamic  stall  vortex.  The  growth  of  the  second  harmonic  indicates  the  flow  is 
highly  nonlinear. 

The  various  turbulence  models  produce  different  sizes  of  dynamic  stall  vortex  and 
separated  flow.  An  examination  of  the  loci  of  the  flow  reversal  point  (xs)  on  the  upper 
surface  of  the  airfoil  (Fig.  49)  shows  that  at  any  given  instant  of  time,  the  extent  of  reversed 
flow  varies  widely.  The  B-L  model  produces  the  smallest  extent  of  reversed  flow  over  most  of 
the  cycle  whereas  the  B-B  model  produces  the  largest  extent  of  reversed  flow.  The  position 
indicating  0°  phase  denotes  the  mean  angle  of  oscillation.  It  is  apparent  from  this  figure 
that  the  B-B  model  completes  the  recovery  process  on  the  upstroke  only  when  it  reaches 
approximately  the  mean  angle.  For  the  large  part  of  the  cycle,  from  15°  upstroke  to  15° 
downstroke,  the  RNG,  B-B,  and  S-A  models  predict  nearly  the  same  extent  of  reversed 
flow.  The  massive  reversed  flow  regions  clearly  increase  the  unsteady-lag  effects. 

Figure  50  shows  a  four-part  figure  of  the  instantaneous  flow  pictures  of  the  extent 
of  the  dynamic  stall  vortex  at  different  times  of  the  oscillatory  cycle  for  the  B-B  model. 
The  flow  which  separates  in  the  trailing-edge  region  during  the  upstroke  cycle  continues  to 
increase  in  extent  by  moving  the  flow  reversal  point,  5,  upsteam  towards  the  leading-edge 
region.  As  seen  in  this  figure,  the  dynamic  stall  vortex  has  peaked  in  in  its  size  around  17° 
downstroke  and  from  then  on  it  prepares  to  shed  by  moving  the  flow  reversal  point  away 
from  the  leading  edge. 

Another  way  of  demonstrating  the  performance  of  the  various  turbulence  models  is 
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through  the  examination  of  instantaneous  flow  pictures  at  any  one  given  phase  in  the  os¬ 
cillatory  cycle.  For  example,  Fig.  51  shows  instantaneous  steamline  particle  flow  pictures 
for  all  turbulence  models  corresponding  to  the  instant  when  the  airfoil  is  at  16°  on  the 
downstroke.  The  three  models,  RNG,  B-B,  and  S-A  have  very  similar  dynamic  stall  vortex 
structures  at  this  instant.  The  B-L  model  has  a  complex  pattern  with  primary  and  sec¬ 
ondary  vortices,  whereas  the  J-K  model  has  already  shed  the  primary  vortex  at  this  time. 
All  models  produce  multiple  vortices  at  slightly  different  times  during  the  downstroke.  The 
B-L  model  also  produces  a  small  bubble  in  the  leading  edge  region,  identified  as  51,  when 
the  airfoil  has  reached  14.5°  in  the  upstroke.  This  bubble  stays  distinctly  separate  from  the 
region  containing  the  dynamic  stall  vortex  and  eventually  merges  with  this  at  about  15° 
downstroke.  The  locus  of  the  reverse  flow  points  shown  in  Fig.  49  is  for  the  point  marked 
5  in  Fig.  51a. 

The  instantaneous  streamline  pictures  of  Fig.  51  are  used  only  for  qualitative  com¬ 
parison  of  different  turbulence  models.  In  fact,  such  a  representation  of  flow  appears  to 
be  misleading  as  it  does  not  depict  the  correct  picture  of  the  unsteady  flow.  Visualization 
of  the  same  flow  using  streaklines  show  a  very  different  flow  picture.  (The  steaklines  are 
computed  using  a  Lagrangian  description  for  the  fluid  particles  through  UFAT75  program.) 
For  example,  Fig.  52  shows  the  same  flow  that  is  pictured  in  Fig.  51d  at  16°  downstroke 
for  the  B-B  model.  As  seen  here,  the  flow  pattern  and  the  details  shown  by  streaklines  are 
phenomenally  different  compared  to  Fig.  51d.  The  large-scale  dynamic-stall  vortex  (VI), 
the  pairing  of  vortices  downstream  of  trailing-edge  (V2),  and  a  diffused  pair  of  vortices 
further  downstream  of  this  (V3)  is  something  that  is  not  apparent  from  the  instantaneous 
streamline  patterns  of  Fig.  51d.  Therefore,  it  is  necessary  to  be  cautious  in  interpreting 
instantaneous  steamline  patterns  of  unsteady  flowfield. 

Another  case  of  deep-stall  considered  is  for  the  mean  angle  of  a0  =  17°.  This  flow 
was  calculated  using  both  10,000  and  20,000  time-steps  per  oscillatory  cycle.  It  appears 
that  for  this  deep-stall  case  at  least  15,000  time-steps  are  needed  to  capture  the  important 
details  of  the  flow.  The  results  of  unsteady  airloads  presented  in  Fig.  53  are  calculated 
using  20,000  time-steps.  Comparison  of  calculations  with  experiments  show  that  all  models 
predict  the  lift  hysteresis  fairly  accurately,  although  the  RNG  and  the  J-K  models  show 
oscillatory  behavior  during  the  downstroke.  The  calculated  results  are  slightly  shifted 
from  the  experimental  data  and  every  model  predicts  separation  at  different  instant  on  the 
upstroke.  Nevertheless,  they  all  reproduce  the  details  of  the  lift  hysteresis  correctly.  The  J- 
K  model  predicts  Cd  and  Cm  loops  that  are  incorrect  both  for  the  upstroke  and  downstroke. 
The  S-A  model  although  has  the  right  trends  for  drag  and  pitching-moment,  it  produces 
separation  too  early,  with  the  result  it  underpredicts  the  peak  drag  and  pitching-moments. 
The  RNG  and  the  B-B  models  calculate  all  the  three  components  fairly  accurately.  The 
RNG  model  shows  oscillatory  behavior  for  the  drag  and  pitching-moment  curves  also.  The 
results  for  all  the  models  are  shifted  compared  to  the  experimental  data.  The  flowfield 
on  the  downstroke  part  of  the  cycle  is  very  complicated.  Although  it  is  shifted  from  the 
experimental  data,  it  appears  that  the  B-B  model  has  demonstrated  superior  performance 
for  this  flow  condition  in  predicting  ail  unsteady  airloads  correctly. 
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It  is  important  to  know  how  much  each  model  costs  computationally.  Figure  54  shows 
the  relative  costs  of  these  turbulence  models  averaged  over  all  the  calculations  performed. 
The  B-L  model  is  used  as  the  reference  as  this  is  the  most  commonly  used  of  all  models 
and  least  expensive.  As  seen,  the  B-B  model  is  the  most  expensive  model  and  costs  about 
2.5  times  the  cost  of  B-L  model.  The  S-A  model  closely  follows  the  B-B  model.  There  is 
no  clear  choice  of  any  one  single  model  that  has  consistently  superior  performance  for  all 
flow  regimes.  The  one  that  comes  close  to  this  from  this  evaluation  is  the  RNG  model.  It 
is  also  cheaper  to  run  and  costs  as  much  as  the  B-L  model.  The  B-L  model  accounts  for 
15.6%  of  the  total  CPU  of  the  numerical  code. 

9.4  Three-Dimensional  Deep-stall  case:  a(f)  =  15°  +  4.2°  sin(Kt) 

The  TURNS  CFD  code  was  modified  for  calculating  the  unsteady  flowfields  of  an  oscil¬ 
lating  wing.  Some  preliminary  calculations  were  done  using  Baldwin-Lomax26  turbulence 
model  for  the  deep-stall  case  of  a(t)  =  15°  +  4.2°  Sin(Kt).  Just  like  in  the  two-dimensional 
case,  the  Baldwin-Lomax  model  was  found  to  be  inadequate  to  resolve  the  separated  flow 
and  the  model  totally  underpredicted  the  extent  of  separation.  Some  findings  of  this  study 
was  presented  at  the  International  Dynamic  Stall  Workshop.  The  inadequacy  of  the  simple 
algebraic  model  to  resolve  the  unsteady  separated  flow  necessitated  a  detailed  study  of 
using  alternative  turbulence  models  that  are  more  accurate.  The  results  of  this  study  is 
described  in  the  previous  sections. 

10.  CONCLUDING  REMARKS 

The  unsteady,  two-dimensional  flowfield  of  an  oscillating  NACA  0015  airfoil  is  calcu¬ 
lated  using  an  implicit,  finite-difference  numerical  method  for  the  solution  of  the  Navier- 
Stokes  equations  with  intent  to  evaluate  the  accuracy  of  five  widely  used  turbulence  models 
to  calculate  the  unsteady  separated  flows  of  dynamic  stall.  Several  unsteady  flow  conditions 
corresponding  to  attached  flow,  light-stall,  and  deep-stall  cases  of  an  oscillating  wing  experi¬ 
ment  were  chosen  as  test  cases  for  calculations.  The  comparison  of  results  with  experiments 
show  that  the  RNG,  the  J-K,  and  the  S-A  models  predict  lift,  drag,  and  pitching-moment 
hysteresis  that  are  in  good  agreement  for  unsteady  attached  flow.  All  three  models  over¬ 
predict  the  extent  of  separation  for  the  light-stall  case  and  therefore  the  airloads  have  good 
agreement  only  for  the  upstroke.  They  have  good  qualitative  agreement  for  lift  and  drag 
hysteresis  for  the  downstroke,  but  they  fail  to  predict  pitching-moments  correctly.  The 
B-L  model  performs  very  poorly  even  for  the  attached  flow.  The  B-B  model  also  performs 
poorly  in  predicting  lift  and  pitching-moment  for  this  attached  flow  case.  For  the  light-stall 
case,  it  predicts  lower  lift  than  experiment  for  the  downstroke  because  of  slow  recovery.  But 
the  drag  and  pitching-moment  are  better  predicted  than  other  models.  For  the  deep-stall 
cases,  the  RNG,  the  B-B,  and  S-A  models  all  predict  qualitatively  correct  airload  hystere¬ 
sis,  but  the  RNG  model  yields  oscillatory  behavior  during  the  downstroke  for  the  extreme 
deep-stall  case.  The  J-K  results  are  not  even  qualitatively  correct  on  the  downstroke  for  the 
deep-stall  cases.  Overall,  the  RNG  model  provides  significant  improvement  over  the  B-L 
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model  with  no  additional  computational  cost.  Among  the  models  considered  here,  the  B-B 
model  is  the  most  expensive  and  costs  about  2.5  times  the  cost  required  to  run  the  code 
with  the  B-L  or  RNG  models.  The  B-L  model  accounts  for  15.6%  of  the  total  CPU  time 
required  to  run  the  numerical  code.  Finally,  visualization  of  unsteady  flowfields  by  means 
of  instantaneous  streamlines  will  give  a  misleading  representation  of  the  flow.  A  Lagrangian 
description  of  particle  motion  using  streaklines  is  perhaps  more  cumbersome  but  better. 
Such  a  description  gives  a  totally  different  flow  picture  compared  to  the  streamline  pattern. 

A  modified  version  of  the  TURNS  code  with  Bald win- Lomax  turbulence  model  was 
used  to  calculate  the  unsteady  flowfields  of  an  oscillating  wing  for  the  deep-stall  case.  Just 
like  in  two-dimensional  case,  the  simple  algebraic  model  was  found  to  be  inadequate  and 
produced  very  little  separation  and  poor  agreement  for  the  results  with  experiments. 
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Table  1 

Calculated  and  experimental  hover  performance  parameter*  tor  the 
UH-60A  Blade  Hawk  and  BERP  rotor*  for  Mttp  =  0.628 


Rotor 

Methods  of  results 

9c 

&r la 

Cq/O 

FM 

UH-60A  model  rotor 

Experiments 

Restricted 

Information 

B9 

mi 

0.73 

N-S  results 

9* 

04)64 

0.0068 

0.73 

N-S  results  with 
new  faritetd  BC 

9* 

04)67 

04)069 

0.76 

UH-60A  equivalent 
rectangular  blade  rotor 

N-S  results 

m 

0.076 

0.69 

UH-60A  linear  twist 
blade  rotor 

N-S  results 

9* 

0.063 

0.0072 

0.68 

UN-60A  model  rotor 

N-S  results 

13* 

0.140 

04)149 

0.72 

BERP  rotor 

N-S  results 

9* 

04)96 

0.0100 

0.72 

N-S  results 

13* 

0.162 

0.0202 

0.76 

Table  2 

Throat  and  power  coefficients  for  the  UH-60A  and  BERP  rotors 


Rotor 

0c 

Or 

Co. 

Cq 

UH-60A 

Black  Hawk 
model  rotor 

9* 

0.007038 

0.0004857 

0.0000867 

04)00572 

13* 

0.011821 

0.0011866 

0.0000694 

0.001256 

BERP  rotor 

9* 

0.010741 

0.0009638 

0.0001291 

0.001093 

13* 

0.017628 

0.002072 

0.0001288 

0.002201 
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Fig.  1  A  C-H  cylindrical  grid  topology  for  a  two-bladed  rotor;  a)  view  in  the  plane 
of  the  rotor,  and  b)  isometric  view  showing  the  grid  boundaries  for  a  single  blade. 
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-  Navier-Stokes  captured  wake  results  -  Present 

- -  Navter-Stokes  prescribed  wake  results  -  Ret.  9 

•  O  Experimental  data  -  Ret.  31 


(b) 


Leading  Edge 


Fig.  5  Computed  surface  particle  flow  detail  highlights  the  shock-induced  bound- 
ary  layer  separation  for  the  flow  conditions  of  a)  MtiP  =  0.877,  6C  —  8°,  and 
Re  =  3.93x10®,  and  b)  Mtip  =  0.794,  6C  =  12°,  and  Re  =  3.55x10®. 


“ “■  **tip  ■  0  877 1  Navier-Stokas 
—  •  —  Mtjp  »  0.44  I  Fin#  Grid  Results 
□  Mtjp  ■  0.877 
O  Mtfp.0.44 

K  Experimental  rang*  0.44  «  Mpp  «  0.877 


Experiments 
(Ret.  31) 


Fig.  6  Comparison  of  bound  circulation  distribution  for  the  case  of  collective  pitch 
6C  =  8°  with  tip  speeds  of  Miip  =  0.44  and  0.877. 
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Blade  Motion 


Fig.  7  Calculated  tip  vortex  particle  flow  details  showing  the  near-field  view  for  the 
flow  condition  of  Mup  =  0.794,  0C  =  12°,  and  Re  =  3.55xl06. 


Blada  motion 


domain 


Fig.  8  Calculated  tip  vortex  trajectory  for  the  flow  conditions  of  A/tip  =  0.44, 
9e  =  8°,  and  Re  =  1.92xl06;  views  show  (a)  contraction  of  wake  (looking  from  the 
top),  (b)  captured  tip  vortex  path  and  its  vertical  descent,  and  (c)  comparison  of 
the  calculated  trajectory  over  450°  of  vortex  age  with  experments. 
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— —  Navwr-Stokas  captured  wake  results  -  Fins  gnd  (217  <71  x  61) 

- Navier-Stokas  captured  wake  results  -  Coarsa  gnd  (109  x  36  *  31) 

•  O  Expanmental  data  -  Rat.  31 


-  Navier-Stokes  captured  wake  results  -  Present 

- Euler  captured  wake  results  -  Present 

•  O  Experimental  data  -  Ref.  31 


=  0.877,  9C  =  8°,  and  Re  =  3.93xl06. 
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(a)  MEM*  rater  MM 


(b>  UH4QA  Nook  Hvak  rater  Meoa 


Fig.  1 1  Planform 
blades. 


views  and  surface  grids  for  the  a)  BERP  and  b)  UH-60A  rotor 


Fig.  12  Radial  twist  distribute 


ns  for  the  UH-60A  and  BERP  rotor  blades. 
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Bottom 

boundary 


Fig  13  View  of  a  C-H  grid  for  a  single  BERP  blade  showing  outer  boundaries  and 
grid  in  the  plane  of  the  blade. 
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-2.8 


-  Nevier-Stokes  results 

•  O  Experiments  (Ref.  6) 


Fig.  14  Surface  pressure  distributions  at  various  span  stations  for  the  UH-60A 
blade  compared  with  experiments.  Mup  =  0.628,  Re  =  2.75  x  10®,  and  calculated 
Ct/o1  —  0.084. 
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Fig.  15  Sectional  thrust  loading  for  the  UH-60A  blade,  equivalent  UH-60A  rectan¬ 
gular  blade,  and  the  UH-60A  blade  with  linear  twist  distribution  in  the  tip  region. 
Map  =  0.628,  0C  =9°,  and  Re  =  2.75  x  106. 


Fig.  16  Comparison  of  calculated  sectional  thrust  loading  for  the  UH-60A  blade  for 
two  thrust  conditions.  Mtip  =  0.628,  and  Re  =  2.75  x  106. 
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Separation 


Collective  pitch  a  13* 


Fig.  17  Upper  surface  particle  flow  (skin  friction)  patterns  of  UH-60A  blade  for 
Bc  =  9°  and  13°  cases  at  Mtip  =  0.628  and  Re  =  2.75  x  106. 


(a)  •  degrees  (b)  13  degrees 


Fig.  18  Close-up  view  of  surface  particle  flow  in  the  tip  region  for  the  UH-6A0  blade 
at  6C  =  9°  and  9C  =  13°.  Mtip  =  0.628,  and  Re  =  2.75  x  106. 
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MACH  CONTOURS  FOR  UH-60  BLADE 

Mtip  =  0.628,  6C  =  13°,  Re  =  2.75  x  106 


View  at  section  A-A 


Fig.  19  Mach  contours  for  the  UH-60A  blade  for  the  conditions  of  Mup  =  0.628, 
Bc  =  13°,  and  Re  =  2.75  x  106.  Views  show  upper  surface  contours  at  L  =  8  (away 
from  the  wall),  and  cross-sectional  plane  at  r/R  =  0.94. 
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kiMi 


CQ,  =  1.15  C^/V2 


Fig.  20  Comparison  of  calculated  and  experimental  hover  performance  for  four- 
bladed  UH-60A  Blackhawk  and  BERP  rotors. 
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-  Ideal 

- Co,  = 

O  UH-60A-CQp|,  -I 

•  UH-60A-Cq 

*  l)H-60A-Eq.  ract.  Made 

▼  UH-SOA-Unear  twist  btada  results 

□  BERP-Cq 

■  BERP-Cq 

UH-60A  modal  experiments  (Ref.  6) 


Fig.  21  Comparison  of  calculated  and  experimental  Figure  of  Merit  for  four-bladed 
UH-60A  Blackhawk  and  BERP  rotors. 
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Fig.  22  Surface  pressure  distributions  for  the  BERP  rotor  blade  at  different  span 
stations;  Muv  =  0.628,  6C  =  13°  and  Re  =  2.94  x  10®. 
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Fig.  23  Radial  circulation  distributions  for  the  UH-60A  and  BERP  blades  at  a 
collective  pitch  setting  of  6C  =  9°  with  Mtiv  =  0.628. 


Fig.  24  Comparison  of  radial  thrust  loadings  for  the  BERP  blade  at  0C  =  9°  and  13° 
with  Mtip  =  0.628  and  Re  =  2.94  x  106. 


Fig.  25  Bound  circulation  distributions  for  the  BERP  blade  for  the  conditions  of 
Fig.  24. 
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Shock-induced 

separation 


Flow  pettern  In  the  notch  region  Pressure  contours  at  A-A 


Fig.  27  Upper  surface  particle  flow  pattern  and  pressure  contours  at  th 
region  for  the  BERP  blade  at  M,ip  =  0.628,  6C  =  13°,  and  Re  =  2.94  x  10' 
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Fig.  29  Schematic  of  new  farfield  boundary  conditions  for  the  hovering  rotor. 


Fig.  30  Sectional  thrust  distributions  for  the  UH-60A  blade  with  old  and  new  farfield 
boundary  conditions  compared  with  experimental  data.  Map  —  0.628,  6C  —  9  »  and 
Re  =  2.75  x  106. 
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Fig.  31  Calculated  wake  contraction  and  descent  for  the  4-bladed  UH-60A  Black 
Hawk  rotor.  Views  show  (a)  planform  view  of  wake  contraction  looking  from  top; 
(b)  wake  contraction  and  descent  looking  from  side;  and  (c)  comparison  calculated 
results  with  experimental  data  for  the  flow  condition  of  Mtip  =  0.628,  9C  =  9°,  and 
Re  =  2.75  x  106. 


Fig.  32  Comparison  of  instantaneous  surface  pressures  for  a  nonliftng  two-bladed 
rotor  in  forward  flight.  Mtip  =  0.6,  p  =  0.2,  and  Re  =  2.22xl06  at  r/R  =  0.893. 


Fig.  33  Instantaneous  surface  pressures  of  a  nonlifting  rotor  in  forward  flight  at 
different  azimuthal  locations  compared  with  experiments.  =  0.8,  n  =  0.2,  and 
Re  =  2.89xl06  at  r/R  =  0.893. 


- 3-D  Navier-Stokes  results 

•  o  Experiments  (Caradonna  et  al.) 


Fig.  34  Instantaneous  surface  pressures  during  transonic  parallel  BVI  at  =  174.5° 
azimuth  compared  with  experiments.  Mtip  =  0.8,  n  =  0.2,  Re  =  2.89xl06,  f  = 
0.177,  x0  =  0.0,  z0  =  -0.4,  at  r/R  =  0.893. 
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Fig.  35  Pressure-time  histories  of  high-speed  impulsive  noise  of  a  hovering  rectan¬ 
gular  blade  rotor  in-plane  at  r/R  =  3.09  and  at  a  blade-tip  Mach  number  of  (a) 
Mtip  =  0.90,  and  (b)  Mtl>  =  0.95. 
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■  Wind  tunnel  (Ref.  14) 
A  Flight  test  (Ref.  14) 


Fig.  36  Pressure- time  histories  of  high-speed  impulsive  noise  for  an  advancing  rotor 
in-plane  at  1.8  rotor  diameters  directly  in  front  of  the  blade  at  Mtip  =  0.896,  p  = 
0.348,  and  Re  =  2.17x10®  compared  to  wind  tunnel  and  flight  test  data. 
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TABLE  3:  QUASI-STEADY  AIRLOADS 


Mnii  angt* 

a.  13* 
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1  — 

Airloads 
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Cd 
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a0208 

aoi4o 

1.624 
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J-K  modal 
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-0.0119 
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1.15 

0.0331 
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0.151 
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Fig.  37  Quasi-steady  surface  pressure  distributions  of  NACA  0015  airfoil  compared 
to  experiments  (Ref.  19)  at  flow  conditions  of  Moo  —  0.29,  Re  =  1.95x10®;  a) 
o  =  11°;  b)  a  =  17°. 


Fig.  38  Comparison  of  calculated  unsteady  airloads  hysteresis  of  oscillat¬ 
ing  wing  for  different  turbulence  models  with  experiments.  Moo  =  0.29, 
a(t)  =  4°  +  4.2 °sin(Kt),  k  =  0.1,  Re  =  1.95xl06. 
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Fig^  39  Comparison  of  harmonic  components  of  unsteady  pressures 
urbulence  models  at  the  expenmental  conditions  of  Fig.  38. 
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Fig.  40  Comparsion  of  unsteady  airloads  hysteresis  for  J-K  model  with  differ¬ 
ent  constant  time-steps  for  the  oscillating  cycle  and  experiments.  Moo  =  0.29, 
a(t)  =  11°  +  4.2° sin(Kt),  k  =  0.1,  and  Re  =  1.95xl06. 
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Fig.  41  Effect  of  numerical  dissipation  on  the  solution  accuracy  with  J 
the  experimental  conditions  of  Fig.  40. 


K  model  for 
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Fig.  42  Effect  of  grid  sensitivity  on  the  solution  accuracy  with  J-K  model  for  the 
experimental  conditions  of  Fig.  40. 
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Fig.  43  Effect  of  grid  sensitivity  on  the  solution  accuracy  with  B-B  model  for  the 
experimental  conditions  of  Fig.  40. 
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Fig.  44  Comparison  of  unsteady  airloads  hysteresis  for  different  turbulence  models 
with  experiments  for  the  flow  conditions  of  Fig.  40. 
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Fig.  45  Comparison  of  harmonic  components  of  unsteady  pressures  for  the  J-K  and 
B-B  turbulence  models  at  the  flow  conditions  of  Fig.  40. 
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flow  conditions  of  Fig.  46  using  the  B-B  turbulence  model. 


Fig.  48  Calculated  harmonic  components  of  unsteady  pressures  for  the  flow  condi 
tions  of  Fig.  46  using  the  B-B  turbulence  model. 
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Fig.  49  Locus  of  flow  reveral  point  during  the  oscillatory  cycle  for  different  turbu 
lence  models  for  the  flow  conditions  of  Fig.  46. 


Fig.  50  Instantaneous  steamline  patterns  at  four  different  time  instants  during  the 
oscillatory  cycle  for  the  B-B  turbulence  model  and  for  the  flow  conditions  of  Fig.  46. 
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Fig.  52  Instantaneous  streakline  pattern  for  flow  at  16°  downstroke  using  the 
Baldwin-Barth  model  and  for  the  flow  conditions  of  Fig.  51. 
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Fig.  53  Comparison  of  unsteady  airloads  hysteresis  for  different  turbulence 
models  with  experiments.  =  0.29,  k  =  0.1,  Re  =  1.95x10®,  and 

a(t)  =  17°  +  4.2°sin(Kt).  98 
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Fig.  54  Comparison  of  relative  costs  of  computations  with  different  turbulence 
models. 
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Euler  Calculations  of  Unsteady  Interaction  of 
Advancing  Rotor  with  a  Line  Vortex 
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Tbc  unsteady,  tbree-dlmeasioaai  flowfield  of  a  helicopter  rotor  Made  ia  forward  flltbt  eocouaiering  a 
coacentrated  line  vortex  is  cakalated  usiat  aa  iaapikit.  finite  difference  aasaerkal  procedure  for  the  solution  of 
Euier  equations,  a  prescribed  vortex  method  is  adopted  to  preserve  the  structure  of  the  interaction  vortex.  The 
test  cases  considered  for  compulation  correspond  to  the  two-bladed  model  rotor  experimental  conditions  of 
Caradonaa  et  al.  and  consist  of  paiaHei  aad  oblique  interactions.  Comparison  of  numerical  results  with  the  test 
data  show  good  agreement  for  the  surface  pressures  for  both  parallel  aad  oblique  interactions  al  subsonic  aad 
traasoak  flow  conditions.  The  results  indicate  that  the  subsoak  parallel  Made-vortex  interaction  is  nearly 
two-dimeasioaai-iike  aad  the  aasteady  time  lag  effects  appear  to  he  negligible.  However,  both  the  three-diama- 
sioaal  aad  unsteady  time-lag  effects  are  found  to  he  important  under  supercritical  flow  conditions  and  these 
effects  are  accentuated  la  the  presence  of  traasoak  shocks.  Ia  contrast,  the  oblique  Made-vortex  interaction  is 
unsteady  and  three  dimensional  al  both  the  subsoak  aad  traasoak  flow  conditions. 


Nomenclature 

a  =  speed  of  sound 

a„  =  vortex  core  radius,  see  Eq.  (8) 

C  =  characteristic  length  scale,  chord  of  the  rotor  blade 

Ck  =  constant,  see  Eq.  (8) 

Cp  =  pressure  coefficient  based  on  local  dynamic 
pressure 

Cw  =  chord  of  the  vortex  generating  wing 
e  =  total  energy  per  unit  volume 

f.G,  ft  ~  flux  vectors 

J  =  Jacobian  of  the  coordinate  transformation 

Mm  =  freest  ream  Mach  number,  forward  speed  of  the 
rotor 

Mtip  -  t>P  Mach  number  of  the  rotor  blade 
p  =  static  pressure 

Q  =  vector  of  conserved  flow  variables 

R  =  rotor  radius 

r  =  radial  distance  from  the  vortex  center 

rB  =  rotor  reference  station  normalized  by  R 

(R(f)  =  rotational  matrix,  see  Eq.  (4) 

U,  V,  W  =  contravariant  velocity  components 
u„  =  forward  flight  speed 

u,  v,  w  =  velocity  components  in  physical  space 

v,  =  tangential  velocity  of  the  vortex 

x,  y,  z.  t  =  inertial  coordinates 

x0,  z„  =  vortex  offset  position  relative  to  the  rotor  axis 
x,.z,  =  distance  from  blade  leading  edge  to  the  line  vortex 

x,  y.  z.  t  =  blade-fixed  coordinates 
7  =  ratio  of  specific  heats 

T  =  vortex  strength 

f  =  dimensionless  strength  of  the  vortex,  T/(amC) 

u  =  advance  ratio,  um/QR 

{.  Vf  f.  r  =  generalized  curvilinear  coordinates 

p  =  density 

♦  =  azimuthal  angle 

Q  =  angular  velocity  of  the  rotor 
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Subscripts 

»  =  freestream 

o  =  interacting  vortex  flow  field 

Introduction 

HE  accurate  simulation  of  the  flowfleld  of  a  helicopter 
rotor  is  still  one  of  the  most  complex  and  challenging 
problems  of  applied  aerodynamics.  This  is  true  in  spite  of  the 
availability  of  the  present-day  supercomputers  of  Cray-2  class 
and  improved  numerical  algorithms.  The  major  reason  for 
this  is  that  the  flowfield  of  a  rotor  in  forward  flight  is  very 
complex.  It  is  highly  three  dimensional,  unsteady,  and  vis¬ 
cous,  with  pockets  of  transonic  flow  near  the  blade  lips  on  the 
advancing  blades  and  regions  of  dynamic  stall  on  the  retreat¬ 
ing  blades.  In  addition,  the  blades  also  shed  complex  vortical 
wakes.  The  concentrated  tip  vortices  shed  by  these  blades 
generally  have  a  close  encounter  with  the  following  blades. 
Such  a  close  encounter  of  the  force-free  concentrated  vortices 
with  the  rotor  blades  is  often  the  cause  of  unsteady  load 
fluctuations  and  impulsive  noise.  These  and  other  complex 
problems  associated  with  such  a  flowfield  are  delineated  in  a 
schematic  of  a  helicopter  in  forward  flight  in  Fig.  1. 

The  spiralling  vortex  sheet  emanating  from  each  of  the 
blades  of  the  rotor  has  a  profound  influence  on  the  perfor¬ 
mance  of  the  helicopter.  It  not  only  alters  the  effective  pitch 
angle  of  the  blades  and  thus  the  airloads,  but  also  produces 
highly  nonlinear  interactions  of  the  vortex  with  the  rotor 
flowfield.  It  is  possible  that  such  interactions  might  produce 
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Fig.  1  Schematic  of  the  complex  flowfieM  of  a  helicopter  ia  forward 

night. 
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vortex-induced  boundary-layer  separation  which  results  in  a 
sudden  loss  of  lift  and  increase  in  drag.  An  accurate  simula¬ 
tion  of  the  rotor  flow  field,  therefore,  must  consider  the  in¬ 
duced  effects  of  the  vortex  wake  including  the  blade-vortex 
interactions  (BV1). 

Numerical  simulations  of  vortex  wakes  are  being  attempted 
only  recently  as  bigger  and  faster  computers  have  become 
available. 1  This  and  other  investigations  have  had  some  lim¬ 
ited  success  to  date.  On  the  other  hand,  carefully  measured 
experimental  data  on  the  three-dimensional  BVI  have  been 
made  available  recently.2-3  But  much  of  the  progress  in  model¬ 
ing  these  blade-vortex  interactions  has  been  hampered  by  the 
lack  of  development  of  theoretical  and/or  numerical  tech¬ 
niques  to  preserve  the  structure  of  the  concentrated  vortices  in 
the  flowfield  without  significant  diffusion.  The  study  of 
blade-vortex  interaction  has  been  the  subject  of  numerous 
recent  research  papers.4'13  These  studies  have  considered  dif¬ 
ferent  methods  of  preserving  the  structure  of  the  concentrated 
vortex  while  convening  in  the  flowfield.  Among  these,  the 
vortex- fitting  (also  called  the  prescribed-vortex  or  perturba¬ 
tion  method  or  split-potential  formulation)  technique  has  been 
demonstrated  to  be  a  very  effective  method  in  preserving  the 
vortex  even  when  the  computational  grid  is  sparse.  The  major 
drawback  of  this  method  is  that  it  does  not  allow  for  the 
distortion  of  the  vortex  core.  For  interactions  where  the  vortex 
does  not  impinge  head-on  onto  the  blade,  this  method  has 
proven  to  have  worked  very  effectively  and  economically7’* 
compared  to  a  more  exact  formulation3  that  allows  for  vortex 
distortion. 

The  purpose  of  this  study  is  to  devise  a  numerical  method 
for  the  solution  of  the  Euler  equations  to  calculate  accurately 
the  unsteady  blade-vortex  interactional  flowfield.  In  particu¬ 
lar,  this  paper  will  focus  on  demonstrating  the  ability  to  calcu¬ 
late  an  interaction  flowfield,  which  is  three  dimensional  and 
unsteady,  for  parallel  and  oblique  blade-vortex  interactions  on 
a  model  helicopter  rotor  tested  in  a  wind  tunnel2-1  at  subsonic 
and  transonic  flow  conditions. 

Governing  Equations  and  Numerical  Scheme 

The  governing  partial  differential  equations  are  the  un¬ 
steady  Euler  equations.  For  generality,  the  equations  are 
transformed  from  the  inertial  Cartesian  reference  frame  (x,  y, 
z,  t)  to  the  arbitrary  curvilinear  coordinate  frame  (£,  ij,  f,  r) 
that  moves  with  the  blade,  while  retaining  strong  conservation 
law  form  to  capture  shock  waves.  The  transformed  equations 
can  be  written  as16 

3r«2  - &)  +  d((f-  f0)  +  3, ((3  -  <3„)  +  3 r(#  -  #„)  =  0  (1) 

where 

"  P  1  I"  PU  ' 

(  P“  PuU  +  %xp 

<5  =  -  pv  ,  P~-  pvU  + 1  ,p 

pw  pwU  + 

.<  J  L  UH-i'P. 

■  pv  i  r  pw  - 

J  puV  +  i\xp  ]  puW  +  ixp 

<3=-  pvV  +  i)yp  ,  #*=-  pvH'  +  fjP  (2) 

pwy  +  i)xp  pwW  +  f,p 

.  VH  +  r>,p  J  L  WH-  ( ,p  . 

Here  is  the  vector  of  the  conserved  flow  variables,  namely, 
the  density  p,  the  three  mass  fluxes  pu,  pv,  and  pw  in  the  three 
coordinate  directions,  and  the  total  energy  per  unit  volume  e. 
Similarly,  (J0  is  a  vector  of  the  conserved  flow  variables  corre¬ 
sponding  to  the  solution  of  the  Euler  equations  for  a  pre¬ 
scribed  line  vortex  aligned  with  the  uniform  frees tream  of 
Mach  number  Mm  in  the  y  direction  and  convening  with  the 
flow  as  shown  in  Fig.  2.  The  vectors  P,  <3,  ft,  and  <3«,  /?„ 


represent  the  appropriate  flux  terms  for  the  two  flows,  respec¬ 
tively.  These  flux  vectors,  as  shown,  are  scaled  by  the  Jaco¬ 
bian  J,  e.g.,  Q  =  J-'Q,  etc.  The  contravariant  velocity  com¬ 
ponents  U,  V,  and  W  are  defined  as 

U  =  £,  +  +  £,v  +  £tw 

V  »  +  I.y  +  «ltw 

W  =  f,  +  iju  +  f„v  +  fjW 

In  the  present  formulation,  £  lies  in  the  chordwise  or 
wraparound  direction,  a  is  in  the  spanwise  direction,  and  f  is 
normal  to  the  blade  surface.  The  terms  £,,  £.,,  £,.  £lt  etc.,  are 
the  coordinate  transformation  metrics.14 

The  velocity  components  u,  v,  w,  and  the  pressure  p  are 
related  to  the  total  energy  per  unit  volume  e  through  the 
equation  of  state  for  a  perfect  gas  by 

P«(y-  1)  [e  -  (p/2Xn2  +  v2  +  w2))  (3) 

This  equation  of  state,  together  with  Eq.  (1),  completes  the 
necessary  equations  that  give  the  entire  flowfield  description. 
These  equations  are  nondimensional.  The  reference  length  and 
velocity  scales  used  are  the  chord  of  the  blade  C  and  the 
frees  tream  speed  of  sound  a.,  respectively.  All  length  scales 
are  normalized  by  C  and  time  by  C/am.  The  primitive  vari¬ 
ables  of  Eq.  (1),  namely,  the  density  p,  the  mass  fluxes  pu,  pv, 
pw,  and  the  energy  per  unit  volume  e,  are  normalized  by  the 
freestream  reference  quantities  made  up  of  p.  and  a.  and  the 
pressure  p  by  ypm. 

In  the  present  formulation,  Eq.  (1)  is  solved  in  the  inertial 
reference  frame  with  the  boundary  conditions  applied  on  the 
routing  blade.  The  terms  u,  v,  and  w  are  the  Cartesian  com- 
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Fig.  2  Schematic  of  an  advaadag  rotor  blade  parting  a  Mm  vertex  in 
a  wind  taaari  experiment  (Kef.  3). 
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ponenis  of  veiocity  in  the  inertial  coordinate  system  u,  y, 
i).  The  inertial  coordinates  ,Y  =  (t.  y,  j.  /)  are  related  to  the 
blade-fixed  coordinates  Xh  -  (x,  y,  z)  through  the  relation 
given  bviT 

•V(x,  v,  c)  =  (R<r)A"„(x.  y.  z),  i  =7  (4) 

where  iR(r)  is  the  rotational  matrix.1*  For  a  rigid  rotor  that 
neglects  the  rotor  disk  attitude,  blade  flapping,  and  pitching 
motions,  the  rotational  matrix  can  be  written  as 

cos  ill  -  sin  Q?  0 

=  -  sinQ?  cosQ?  0  (j) 

J  0  0  1 

L  J 

Here  0  is  the  angular  velocity  of  the  rotor  and  ill  represents 
the  azimuthal  sweep  of  the  rotor  blade. 

The  governing  equations  are  solved  using  a  two-factor, 
implicit.  Unite  difference  numerical  scheme.  The  numerical 
algorithm  that  was  originally  developed  by  Ying  et  al.”  has 
been  modified  to  accommodate  a  prescribed- vortex  perturba¬ 
tion  formulation  for  a  rotating  blade  environment.  The  nu¬ 
merical  scheme  uses  spatial  central  differencing  in  the  it  and  f 
directions  and  upwind  differencing  in  the  f  direction.  The 
two-factored  scheme  containing  the  vortex-fitting  terms  is 
given  by 

[/  +  h6*(A  *)"  +  -  £>,  If] 

x  [/  +  htfjA-r  +  h6j"-  D,  -  A&") 

=  -  a*  i  sfu/? ♦  r  -  (f:  ri  + «{[( f-r-  (F;  n 

+  6,(<S"  -  (5g)  +  -  ft!)] 

-(D,l,  +  £>,  I  ,)(£*-&)  (6) 

where  h  =  At  for  first-order  time  accuracy;  S  is  typically  a 
three-point,  second-order-accurate,  central  difference  opera¬ 
tor;  and  the  operators  6*  and  t/(  are  backward  and  forward 
three-point  difference  operators.  The  time  index  is  denoted  by 
n  such  that  r  =  (/iAr),  and  A<5"  =  0"  * 1  -  Q"-  The  flux  vector 
P  has  been  split  into  P*  and  P\,  according  to  its  eigenval¬ 
ues,20  and  the  Jacobian  matrices  A  *,  B,  and  C  result  from  the 
local  linearization  of  fluxes  about  the  previous  time  level.16  In 
writing  Eq.  (6),  it  is  assumed  that  A*  *  A  *,  6„  «■  6 ,  and 
C„  =»  C  where  the  Jacobian  matrices  A* ,  B0,  and  C0  corre¬ 
spond  to  the  prescribed-vortex  flowfield.  In  the  absence  of 
vortex  interaction,  the  prescribed-disturbance  flowfield  re¬ 
duces  to  a  freestream.  The  finite  difference  scheme  described 
in  Eq.  (6)  uses  flux  splitting  in  £  direction  and  central  differ¬ 
encing  in  the  n  and  f  directions.  As  a  consequence,  numerical 
dissipation  terms  A  and  D,  are  used  in  the  *  and  f  directions, 
and  are  given  as  combinations  of  second  and  fourth  differ¬ 
ences.  For  example,  these  terms  in  it  direction  are  given  by 

D,  I,  =  (A/)/-'(l»)lx  +  1*1,  +  l*lt) 

<7*» 

A  I,  =  (A/)/-,(l*lx  +  l*lr  +  lv| I*) 

where  <2  and  <<  are  constants  and  5  is  a  midpoint  operator.  In 
the  vicinity  of  'hock  waves,  the  fourth-difference  terms  can 
cause  oscillations,  so  it  is  desirable  to  drop  these  terms  and 
rely  only  on  the  second-difference  terms. 

The  factored  operators  are  solved  by  sweeping  in  the  ( 
direction  and  inverting  tridiagonai  matrices  with  5x5  blocks 
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for  the  other  two  directions.  Currently,  a  significant  part  of 
the  computational  time  is  taken  to  form  the  plus  and  minus 
Jacobian  matrices  for  the  flux  vector  F  with  this  numerical 
scheme.  However,  this  effort  has  been  reduced  by  computing 
A  *  and  A  at  every  other  point  (in  both  it  and  f  directions) 
and  averaging  to  obtain  the  matrices  at  the  intermediate 
points.  The  numerical  algorithm  is  second-order  accurate  in 
space  and  first  order  in  time  and  the  code  is  vectorized  for  the 
Cray-2  supercomputer. 

A  body-conforming  finite  difference  gnd  has  been  used  for 
the  rectangular  blade  having  an  aspect  ratio  of  7  and  consists 
of  a  warped  spherical  0-0  gnd  topology.  The  flowfield  gnd  is 
generated  using  the  three-dimensional  hyperbolic  grid  solver 
of  Steger  and  Chaussee21  with  proper  clustering  in  the  leading- 
and  traiiing-edge  regions  and  in  the  tip  region.  The  grid  is 
nearly  orthogonal  at  the  surface,  the  spacing  in  the  normal 
direction  at  the  surface  is  chosen  to  be  0.02C,  and  ail  of  the 
calculations  were  done  on  a  21  x  101  x  15  points  grid.  The 
grid  boundary  is  located  at  8  chords  in  all  directions. 

The  boundary  conditions,  both  surface  and  far  field,  are 
applied  explicitly.  The  slip  conditions  use  an  extrapolation  of 
contravariant  velocities  to  the  surface.  The  density  at  the  wall 
is  determined  by  the  zeroth-order  extrapolation.  The  pressure 
along  the  body  surface  is  calculated  from  the  normal  momen¬ 
tum  relation  (see,  e.g..  Ref.  16).  Having  calculated  the  density 
and  pressure,  the  total  energy  is  determined  from  the  equation 
of  state. 

At  the  far-fleld  boundary,  the  flow  quantities  are  either 
fixed  or  extrapolated  from  the  interior  depending  on  whether 
the  flow  is  subsonic  or  supersonic  and  if  it  is  of  inflow-  or 
outflow-type  at  the  boundary.22-21  The  characteristic  velocities 
of  the  Euler  equations  determine  the  number  of  flow  proper¬ 
ties  to  be  specified  to  control  the  reflections  of  waves  from  the 
boundaries.  For  the  subsonic-inflow  boundary,  four  quanti¬ 
ties  are  specified  and  one  quantity  is  determined.  The  four 
specified  here  are  a  Riemann  invariant,  the  entropy,  and  two 
tangential  velocities;  the  quantity  that  is  calculated  is  also  a 
Riemann  invariant.  For  the  supersonic  inflow,  all  flow  quanti¬ 
ties  are  specified.  At  the  subsonic-outflow  boundaries,  only 
one  quantity  is  specified.  For  the  supersonic-outflow  condi¬ 
tion  ail  flow  quantities  are  extrapolated  from  the  interior.  The 
plane  containing  the  blade  root  is  chosen  very  close  to  the 
rotation  axis  of  the  blade  (at  R  =  1 .0 C)  and  is  also  treated  as 
a  far-field  boundary. 

The  concentrated  line  vortex  generated  upstream  of  the 
rotor  by  a  rectangular  wing  is  fixed  in  the  inertial  space  along 
a  line  of  constant  x.  It  is  assumed  to  have  an  analytical 
representation  for  the  cylindrical  velocity  distribution  given  by 


The  constant  C,  in  the  preceding  expression  has  been  deter¬ 
mined  to  be  equal  to  0.8  by  matching  the  peak  tangential 
veiocity  with  the  experimentally  measured  value  at  the  mea¬ 
sured  radial  distance.*  The  constant  a0  is  approximately  equal 
to  the  radius  of  the  viscous  core  of  the  vortex,  and  this  is  equal 
to  the  experimentally  measured  value  of  0.167C.  The  induced 
velocity  field  due  to  this  line  vortex  is  calculated  using  the 
Biot-Savart  law  and  the  pressure  field  is  calculated  by  solving 
the  radial  momentum  equation. 

Results  and  Discussion 

All  of  the  calculations  performed  in  this  study  are  done  in  a 
time-accurate  manner.  Three  test  cases  have  been  chosen  for 
calculations  from  among  the  many  test  conditions  of  the  two 
sets  of  experiments  of  Caradonna  et  al.2J  for  parallel  and 
oblique  blade-vortex  interactions.  Since  one  of  the  purposes  of 
the  experiments  was  to  collect  data  to  validate  numerical 
methods,  the  experimental  apparatus  was  kept  simple  to  ease 
the  representation  by  numerical  methods.  The  rotor  geometry 
consists  of  a  two-bladed  rigid  rotor  of  approximately  14-chord 
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diameter.  The  blades  have  a  rectangular  planform  and  are 
made  up  of  NACA  0012  airfoils  with  no  twist  or  taper.  The 
rotor  blades  are  set  at  zero  collective  pitch  and  are,  therefore, 
essentially  nonlifting  in  the  absence  of  the  vortex  interaction. 
The  interacting  vortex  was  generated  upstream  of  this  rotor  by 
a  lifting  NACA  0015  rectangular  wing.  Parallel  blade-vortex 
interaction  flowfields  were  generated  by  positioning  the  line 
vortex  along  the  y  axis  ( x0  -  0.0).  Similarly,  the  oblique  BVI 
flowfields  were  generated  by  moving  the  position  of  the  line 
vortex  such  that  it  is  still  aligned  along  the  y  axis  with  x0  <0.0. 
This  enables  the  advancing  rotor  blade  to  encounter  the  line 
vortex  continuously  in  the  first  and  second  quadrants  of  the 
azimuthal  travel  as  shown  in  Fig.  2. 

Although  the  main  focus  of  the  present  investigation  is  to 
calculate  the  parallel  and  oblique  blade-vortex  interactions,  it 
is  necessary  to  first  calculate  the  baseline  flowfields  of  rotor- 
alone  configuration  (in  the  absence  of  vortex  interaction)  at 
the  same  freestretun  conditions.  Two  sets  of  results  for  parallel 
blade-vortex  interactions,  corresponding  to  =  0.6  and  0.8 
and  n  =  0.2,  and  one  set  of  results  for  oblique  blade-vortex 
interaction  corresponding  to  M, jp  =  0.763  and  n  =  0.197,  will 
be  presented  in  the  following  sections  and  compared  with  the 
experimental  data. 

Rolor-AioM  Case 

The  calculation  of  the  rotor-alone  flowfield  solutions  serves 
two  purposes.  First,  it  enables  an  understanding  of  the  impor¬ 
tance  of  unsteady  time-lag  effects  in  shock-wave  growth  and 
decay  as  well  as  the  three  dimensionality  of  the  flowfield  of  the 
advancing  rotor  under  transonic  flow  conditions.  Second,  it 
provides  the  baseline  solution  for  starting  the  unsteady  calcu¬ 
lations  of  the  vortex-interaction  flowfield. 

Figure  3  shows  the  instantaneous  surface  pressure  distribu¬ 
tions  for  various  azimuthal  positions  of  advancing  rotor  for 
the  flow  condition  of  M,w  =  0.6  and  n  »  0.2.  Examination  of 
these  results  suggests  that  at  this  reference  station  of  r$  =  0.893 
and  for  the  subcritical  flow  condition,  the  flow  behaves  as  if  it 
is  quasisteady  and  quasi-two  dimensional.  In  spite  of  the 
gradual  increase  of  the  blade-element  relative  speed  in  the  first 
quadrant  and  gradual  decrease  of  the  same  in  the  second 
quadrant,  the  flowfield  appears  to  remain  nearly  the  same  at 
all  azimuthal  locations,  indicating  that  the  unsteady  time-lag 
effects  are  virtually  absent  for  this  flow  condition  and  that  the 
flow  behaves  as  if  it  is  quasisteady.  In  fact,  also  shown  in  Fig. 
3  is  the  quasisteady  surface  pressures  for  one  azimuth  location 
of  *-90  deg,  which  is  in  perfect  agreement  with  the  unsteady 
surface  pressures  at  different  azimuthal  positions.  Also,  the 
comparison  of  these  three-dimensional  results  with  the  two-di¬ 
mensional  results  of  Ref.  8  further  confirms  that  the  two 


Fig.  3  Coagaftma  of  wkihlri  surface  prusani  for  aa  advaadag 
rotor  with  the  sssdsnsSy  surface  prieearn  (Map  -  0.6,  p  -  g.2,  sat 
rg  m  0.103). 


results  are  nearly  identical,  suggesting  that  the  flow  also  be¬ 
haves  as  if  it  is  quasi-two  dimensional. 

At  the  supercritical  tip  flow  conditions,  corresponding  to 
At#,  =  0.8  and  n  =  0.2,  the  basic  rotor  flowfield  is  dominated 
by  the  presence  of  a  strong  shock  wave  on  the  advancing  blade 
over  large  pans  of  the  first  and  second  quadrants.  The  instan¬ 
taneous  surface  pressures  for  the  advancing  rotor  are  shown  in 
Fig.  4  at  the  radial  station  rt  =  0.893  for  different  azimuthal 
positions  of  the  rotor  blade.  The  results  show  good  agreement 
with  the  experiments  for  all  azimuthal  locations.  The  absence 
of  viscosity  and  therefore  the  boundary  layer  causes  the  Euler 
results  to  overpredict  the  shock  strength  and  position  in  the 
first  and  second  quadrants  of  azimuthal  travel.  The  experi¬ 
mental  data  shows  that  the  strong  shock  wave  which  was 
present  at  the  beginning  of  the  second  quadrant  nearly  decays 
by  the  time  the  blade  reaches  *  -  180-deg  position.  In  con¬ 
trast,  the  numerical  method  predicts  that  the  shock  wave 
continues  to  persist  at  least  until  after  it  passes  the  *  -  180- 
deg  azimuthal  location,  in  fact,  it  is  completely  absent  by  the 
time  the  rotor  blade  reaches  *  «  184-deg  position.  Previous 
two-dimensional  Navier-Stokes  calculations  of  the  same  flow* 
had  shown  strong  three  dimensionality  and  unsteady  time  lags 
in  shock-wave  growth  and  decay.  The  two-dimensional  as¬ 
sumption  for  this  flow  essentially  overpredicted  the  shock- 
wave  position  and  strength,  unlike  the  subcritical  flow  condi¬ 
tion.  In  contrast,  the  three-dimensional  Euler  results  seem  to 
follow  the  experimental  observation  correctly. 

Thus  at  this  transonic  flow  condition,  the  flowfield  is  highly 
three  dimensional,  and  exhibits  strong  unsteady  time-lag  be¬ 
havior  in  the  shock-wave  formation  and  eventual  demise.  The 
unsteady  time-lag  characteristic  of  the  flowfield  is  demon¬ 
strated  in  the  surface  pressure  distributions  of  Fig.  4  shown 
for  *  -  60-deg  and  120-deg  azimuthal  positions,  which  shows 
that  the  two  results  are  very  different  from  each  other. 

In  principle,  the  shock  wave  should  attain  its  maximum 
strength  when  the  relative  flow  speed  reaches  a  maximum 
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Fig.  5  Instantaneous  surface  pressure  distributions  at  tbe  reference 
blade  section  for  retreating  blade  (Mug  «  0.8.  a  *  0.2.  and  rg  »  0.893). 


Fig.  6  Comparison  of  quasisteady  and  unsteady  surface  pressures 
for  advancing  rotor  (Mgg  =»  0.8,  g  *  0.2,  sod  rg  “  0.893). 


value  for  ♦  =  90-deg  blade  azimuthal  position.  However,  the 
shock  wave  keeps  growing  in  strength  even  after  the  blade  has 
passed  the  'f'  =  90-deg  azimuthal  position  and  a  weak  shock 
wave  seems  to  persist  even  at  the  ♦  *  180-deg  azimuth.  The 
strengthening  of  the  shock  wave  in  the  first  quadrant  and 
beyond,  and  the  slow  demise  of  it  in  the  second  quadrant, 
suggests  the  existence  of  a  strong  unsteady  time-lag  effect  at 
this  flow  condition.  Such  a  behavior  was  totally  absent  at  the 
subcritical  flow  condition.  The  presence  of  the  shock  wave 
seems  to  introduce  not  only  the  time  lag  in  the  adjustment  of 
the  flow  as  the  blade  sweeps  in  azimuth,  but  also  makes  the 
flow  highly  three  dimensional. 

However,  this  behavior  is  confined  to  the  advancing  side  of 
the  rotor  blade.  In  contrast,  the  flowfleld  of  the  retreating  side 
appears  really  benign  for  this  nonlifting  rotor.  The  lingering 
effects  of  the  shock  wave,  persisting  at  ♦  =  180  deg,  soon  die 
out  as  the  blade  sweeps  into  the  third  quadrant.  Shown  in  Fig. 
5  are  the  surface  pressure  distributions  for  several  azimuthal 
locations  of  the  blade  in  the  third  and  fourth  quadrants  on  the 
retreating  cycle.  Since  the  flowfield  is  basically  subcritical  on 
the  retreating  side,  the  instantaneous  surface  pressures  of  Fig. 
5  appear  to  have  the  same  behavior  as  those  of  Fig.  3. 

The  decay  of  the  strong  time-lag  effects  on  the  retreating 
blade  allows  a  simplification  in  starting  the  unsteady  calcula¬ 
tions,  but  only  if  the  strongest  transonic  region  is  avoided 
initially.  To  illustrate  this  point.  Fig.  6  shows  the  quasisteady 
and  unsteady  surface  pressures  for  the  case  of  Mtip  =  0.8  and 
g  =  0.2  at  the  azimuthal  locations  of  *  =  60  deg  and  120  deg. 
For  this  flow  condition,  the  results  indicate  that  the  quasi¬ 
steady  assumption  overpredicts  the  shock  strength  at  ♦  =  60 
deg  and  underpredicts  at  ♦  =  1 20  deg.  However,  at  ♦  -  0.0-deg 


position,  the  quasisteady  and  unsteady  results  are  nearly  iden¬ 
tical.  even  for  this  flow  condition.  Therefore,  it  is  still  reason¬ 
able  to  start  the  unsteady  calculations  from  the  quasisteadv 
solution  calculated  at  ♦  =  0.0-deg  azimuthal  location. 

Parallel  Blade- Vortex  loteractioa 

During  the  unsteady  three-dimensional  close  encounter  of  a 
curved  tip  vortex  with  a  rotating  blade,  the  helicopter  rotor 
undergoes  a  variety  of  blade-vortex  interactions  depending  on 
the  interaction  angle  between  the  leading  edge  of  the  rotating 
blade  and  the  curved  line  vortex.  These  interactions  are  gener¬ 
ally  unsteady  and  three  dimensional.  One  limiting  case  of  such 
an  encounter,  for  zero  interaction  angle,  is  termed  parallel 
interaction  (see,  c.g..  Ref.  7).  In  the  experimental  configuration 
considered  here,  this  interaction  occurs  around  *  =  180-deg 
azimuthal  position.  For  an  observer  riding  with  the  blade  (at  a 
given  reference  station  along  the  span),  it  appears  as  though 
the  observer  is  passing  a  fixed  vortex  in  the  flow.  For  this 
reason,  this  interaction  is  sometimes  approximated  as  two-di¬ 
mensional  and  unsteady. 

To  calculate  accurately  the  blade-vortex  interaction  flow- 
field,  it  is  necessary  to  preserve  the  vortex  structure  without 
numerical  diffusion.  As  mentioned  before,  one  method  that 
has  been  demonstrated  to  work  effectively  and  economically 
in  achieving  this  is  the  prescribed-disturbance  scheme.7'1  The 
effectiveness  of  this  scheme  is  illustrated  in  Fig.  7,  showing  the 
variation  of  lift  coefficient  as  a  function  of  vortex  position 
during  a  two-dimensional  airfoil-vortex  interaction  calculated 
using  the  Navier-Stokes,  Euler,  and  transonic  small-distur¬ 
bance  methods.  Also  shown  in  this  figure  is  the  Euler  lift 
distribution  calculated  using  a  conventional  (vertex  capturing) 
technique.24  Although  the  two  Euler  solutions  are  computed 
on  the  same  grid,  the  numerical  dissipation  associated  with  the 
finite  difference  grids  progressively  weakens  the  gradients  and 
reduces  the  effective  vortex  strength  in  the  conventional 
method.  Also,  this  numerical  error  is  grid  dependent  (the  finer 
the  grid  the  lesser  the  error);  however,  it  is  completely  absent 
in  the  prescribed-vortex  solutions,  which  are  essentially  grid 
independent.  Hence,  the  prescribed-disturbance  method  has 
been  used  here  for  preserving  the  vortex  structure. 

Subcritical  Interaction 

The  results  of  a  subcritical  parallel  BVI  are  discussed  here. 
This  case  corresponds  to  the  flow  conditions  of  A/,ip  =  0.6, 
p  =  0.2,  and  C*f  =0.133  at  a  blade  reference  station 
rB  =  0.893.  The  interacting  vortex  is  located  at  x„  =  0.0  and 
Zo  =  0.4  along  the  y  axis.  To  calculate  the  BVI  flowfield,  the 
interacting  vortex  is  introduced  into  the  baseline  rotor  solution 
at  the  azimuthal  position  of  ♦  =  120  deg.  The  unsteady  flow- 
field  is  monitored  as  the  blade  is  set  in  motion  to  advance  in 
azimuth.  The  blade-vortex  interaction  effects  peak  around 
♦  =  180-deg  azimuthal  position.  The  instantaneous  surface 
pressure  distributions  are  shown  in  Fig.  8  for  several  az- 


—  — - Euler-nonperturtootton 


Fig.  7  Demons! ratios  of  prescribed-vortex  method  for  preserving 
vortex  structure  in  two-dtoeurioual  elrfoB-vortex  interaction:  NACA 
MAOOf  alrfoN;  Mm  -  0.85.  r  -  0.2.  xad  z,  -  -  0.20. 
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imuthai  locations.  The  present  results  are  compared  with  the 
experimental  data*  and  the  results  from  published  two-dimen¬ 
sional  Navier-Stokes  calculations.* 

The  present  results  are  in  very  good  agreement  with  both  the 
experimental  data  and  the  two-dimensional  Navier-Stokes  re¬ 
sults.  As  seen  from  Fig.  8a  for  ♦  »  178. IS  deg,  the  lift  on  the 
blade,  which  is  initially  zero  in  the  absence  of  vortex  interac¬ 
tion,  is  negative.  Because  of  the  sense  of  its  rotation,  the 
approaching  vortex  introduces  a  downwash  in  the  flow  imme¬ 
diately  ahead  of  the  blade.  The  lift  rapidly  becomes  positive  as 
the  blade  passes  the  line  vortex.  This  crossover  of  lift,  from 
negative  to  positive,  seems  to  occur  when  the  vortex  is  approx¬ 
imately  aligned  with  the  quarter-chord  line  of  the  blade  (corre¬ 
sponding  to  a,  =  0.0).  As  mentioned  before,  the  two-dimen¬ 
sional  approximation  of  this  three-dimensional,  unsteady 
interaction  is,  in  fact,  a  very  good  assumption  for  this  subcrit- 
ical  flow  condition.  The  close  agreement  of  the  viscous  and 
in  viscid  results  suggests  that  the  viscous  effects  are  unimpor¬ 
tant  for  this  relatively  weak  interaction. 

Sapcrcritkal  Interaction 

This  case  corresponds  to  the  flow  condition  of  Mtip  =  0.8, 
/i  >0.2,  and  C*f  =  0.177  at  a  blade  reference  station  of 
rt  >  0.893.  The  interacting  vortex  location  is  the  same  as  in 
the  subcritical  case.  As  in  the  rotor-alone  case  this  flow  condi¬ 
tion.  in  contrast  to  the  subcritical  flow  condition,  exhibits 
strong  unsteady  time-lag  effects  in  shock-wave  growth  and 
decay  for  the  advancing  rotor.  This  feature  seems  to  accentu¬ 
ate  the  three-dimensional  nature  of  the  flow,  and  a  two-di¬ 
mensional  approximation  of  this  flow  overpredicts  the  shock- 
wave  strength  and  location,  as  demonstrated  in  Ref.  8. 

As  before,  to  calculate  the  BVI  flowfieid,  the  interacting 
vortex  is  introduced  at  ♦  -  120-deg  azimuthal  position,  and 
the  time  evolution  of  the  flow  is  monitored  as  the  blade 
advances  in  azimuth.  Figure  9  shows  the  interaction  flowfieid 
results  in  terms  of  instantaneous  surface  pressures  for  several 
azimuthal  positions  of  the  blade.  Also  shown  in  these  figures 
are  the  experimental  data.2  The  present  results  are  in  good 
agreement  with  the  experimental  data  for  all  azimuthal  loca¬ 
tions  of  the  blade,  and  they  capture  the  essential  details  of  the 
flow  sharply,  including  the  shock  waves. 


M  Made  vortex  bwwscdoa  for  Irian  air  caw  (Mw  -•!*!"*  >  IJ, 
C*f  -  *.177,  x.  -  •.§,  *,  -  0.4,  aad  r,  .  0.003). 


Fig.  g  Cowpartwa  of  lartiaiianai  surface  pitmans  daring  prd- 
M  Made  vortex  Interaction  for  takooate  caw  (Af*,  >  M,  p  >  0.2, 
Of  •  0.133,  x,  -  0.0,  U  -  0.4,  aad  ra  -  0  JOS). 


The  effect  of  vortex  interaction  is  to  induce  time-dependent 
aerodynamic  forces  on  the  rotor.  For  example,  as  seen  from 
the  surface  pressure  plots  of  Fig.  9,  the  lift  on  the  blade,  which 
is  initially  zero,  becomes  negative  due  to  induced  downwash 
and  rapidly  increases  to  a  positive  value  as  the  blade  passes  the 
stationary  vortex.  The  peak  effects  of  the  interaction  appear 
to  occur  when  the  blade  leading  edge  is  approximately  above 
the  line  vortex. 

OMIgae  Wade  Vortex  Interaction 

An  interaction  is  called  oblique  when  the  line  vortex  inter¬ 
acts  over  a  large  part  of  the  blade  in  the  radial  direction  at  any 
given  instant  of  its  azimuthal  travel.  In  other  words,  an  ob¬ 
server  fixed  to  the  blade  and  moving  with  it  sees  at  any  given 
instant  different  parts  of  the  line  vortex  interacting  with  the 
blade  as  it  sweeps  from  0  to  180  deg  in  azimuth  (in  the  present 
case).  As  seen  in  the  schematic  of  Fig.  2,  an  oblique  interac¬ 
tion  is  therefore  possible  whenever  x0  is  nonzero.  This  is  in 
contrast  to  the  parallel  blade-vortex  interaction  that  occurs 
when  x„  >  0.0  where  an  observer  located  at  one  radial  station 
sees  the  blade  passing  a  point  vortex  that  is  stationary.  The 
observer  is  therefore  concerned  about  the  influence  of  the 
entire  line  vortex  at  this  particular  radial  station.  Therefore, 
for  an  observer  stationed  on  the  blade,  the  interaction  experi¬ 
enced  of  parallel  and  oblique  BVI  is  uniquely  different. 

The  oblique  interaction  considered  here  corresponds  to  a 
freestream  condition  of  >  0.763  and  p  >  0.197.  The  inter¬ 
acting  vortex  located  at  x.-  -2.13  and  z,  >0.25  has  a 
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Fig.  10  Comparison  of  Instantaneous  surface  pressures  daring 
oblique  blade-vortex  interaction  for  transonic  case  (Afu,  -  0.763, 
a  -  0.1*7,  C*f  -  0.179,  x.  »  -  2.13,  z.  -0.25,  and  r&  -  0.946). 


strength  of  Ckt  =  0.179  for  the  reference  blade  section  at 
rB  =  0.946.  In  contrast  to  the  parallel  interaction,  the  oblique 
interaction  occurs  over  a  larger  azimuthal  sweep  of  the  rotor 
blade,  starting  in  the  first  quadrant  and  completing  in  the 
second  quadrant.  The  peak  interaction  occurs  at  around 
♦  =  20  deg  in  the  first  quadrant  and  around  ♦  =  160  deg  in 
the  second  quadrant.  Therefore,  the  vortex  is  introduced  into 
the  flow  at  ♦  =  0.0  deg  in  this  case,  and  the  evolution  of  the 
unsteady  interacting  flow  is  monitored  for  the  advancing  rotor. 

Figure  10  shows  the  calculated  instantaneous  surface  pres¬ 
sure  distributions  at  several  azimuthal  positions  of  the  blade 
and  the  experimental  data.}  The  comparison  of  results  shows 
overall  agreement  of  surface  pressures.  This  set  of  experimen¬ 
tal  data  appears  to  have  large  scatter.  Even  the  surface  pres¬ 
sures  of  rotor-alone  case  have  similar  scatter.23  The  interaction 
in  the  first  quadrant  occurs  at  subcritical  flow  and  the  agree¬ 
ment  with  data  is  good  at  ♦  =  50  deg.  The  interaction  in  the 
second  quadrant  occurs  under  supercritical  flow  condition. 
The  maximum  relative  blade  tip  speed  in  this  case  (for  ♦  =  90- 
deg  position)  is  0.91  compared  to  a  value  of  0.96  correspond¬ 
ing  to  the  transonic  parallel  interaction  case  and  is  thus  a 
weaker  interaction. 

From  the  surface  pressures  of  Fig.  10  it  can  be  seen  that  the 
lift  on  rotor  blade,  which  is  initially  zero,  becomes  negative 
due  to  the  vortex  influence  (because  of  the  sense  of  rotation), 
increases  to  a  maximum  (negative)  value  in  the  second  quad¬ 
rant,  and  changes  over  to  a  positive  value  at  about  ♦  =  165 -deg 
azimuthal  position  of  the  blade.  The  peak  interaction  appears 
to  occur  around  the  ♦  -  160-deg  azimuthal  position. 


Conclusions 

A  numerical  procedure  is  presented  to  calculate  the  un¬ 
steady,  inviscid.  three-dimensional  flowfieids  of  a  helicopter 
rotor  in  forward  flight  encountering  parallel  and  oblique 
blade-vortex  interactions  in  subsonic  and  transonic  flow  con¬ 
ditions.  Important  flow  features  such  as  unsteady  time-lag 
effects  in  shock-wave  growth  and  demise,  as  well  as  the  impor¬ 
tance  of  three-dimensional  effects,  are  discussed.  Although  it 
is  possible  under  certain  flow  conditions  to  approximate  the 
parallel  blade-vortex  interaction  as  two  dimensional  and  un¬ 
steady,  the  oblique  blade-vortex  interaction,  on  the  other 
hand,  is  strongly  three  dimensional  and  unsteady,  and  there¬ 
fore  cannot  be  approximated  as  two-dimensional  interaction. 

The  numerical  results  are  compared  with  two  sets  of  experi¬ 
mental  data  generated  by  Caradonna  et  al.21  on  a  model 
two-bladed  rotor  in  a  wind  tunnel.  The  present  Euler  results 
show  good  agreement  with  experiments  for  both  the  parallel 
and  oblique  interactions  at  subsonic  and  transonic  flow  condi¬ 
tions.  Thus,  the  numerical  methodology  presented  here  has 
demonstrated  the  ability  to  accurately  calculate  the  flowfield 
of  an  advancing  helicopter  rotor  including  the  effects  of  vor¬ 
tex  interaction. 
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